c     Two cavity pulse compressor using envelope equations.
c     *****************************************************
c     Last modified August 2017

c     This code is provided as a reference document to the PRAB publication
c     "Control and performance improvements of an X-band SLED1 pulse compressor
c      in use for testing accelerating structures at high power."
c     by 
c     Benjamin Woolley and Igor Syratchev
c     CERN, European Organization for Nuclear Research, 1211 Geneva, Switzerland
c     and 
c     Amos Dexter
c     Lancaster University, Lancaster, LA1 4YW, UK

c     The program solves envelope equations for two pulse compressor cavities,
c     combines power with a perfect 3dB Hybrid and allows a phase program
c     on the input as determined by a controller which seeks constant amplitude on
c     the compressed output pulse.

c     To complile use ANSI FORTRAN77 
c     To run the datafile "indata_pha.txt" must be in a place where the program
c     can find it

c     The output files are configured to provide data as required for the
c     published figures.

c23456789012345678901234567890123456789012345678901234567890123456789012

      integer    jpulse, jplot
      parameter (jpulse=8000)
      parameter (jplot=10000)

      real*8 QLA, QEA, QLB, QEB, Q0A, Q0B
      real*8 fcA, fcB, fd, Vd0, fbandwidth
      real*8 pulse_time_mod, flip_time_klys, pulse_time_klys
      real*8 phase_flip, errtime, fac1, fac2

      common/com_imp/QLA, QEA, QLB, QEB, Q0A, Q0B,fcA, fcB, fd, Vd0,
     +               fbandwidth, pulse_time_mod, flip_time_klys,
     +               pulse_time_klys, phase_flip,
     +               errtime, fac1, fac2 

      real*8 hdt, dt6, dt
      real*8 hpi, pi, tpi, recip_root2
      real*8 Vd3, Vd, dV, v1, v2
      real*8 wcA, wcB, wd
      real*8 t0, t2, tstart, ts, time
      real*8 psi0, psi1, psi2, psi, Rpsi, dp1, dp2
      real*8 Rpsi_deg, Psi_deg
      real*8 g1, g2, drr, dri
      real*8 RAx, RAy, RBx, RBy
      real*8 Vx, Vy, Vx0, Vy0
      real*8 DVx, DVy, VAx, VAy, DVAx, DVAy, VBx, VBy, DVBx, DVBy
      real*8 g3, Rx, Ry, Ramp
      real*8 sum, diff, error, error1
      real*8 power_top, err

      real*8 amp1(jplot), amp2(jplot), amp3(jplot), amp4(jplot)
      real*8 tim(jplot)
      real*8 ph1(jplot), ph2(jplot), ph3(jplot), ph4(jplot)
      real*8 ph(jpulse)
      real*8 amp(jplot), ampx(jpulse), dpe(jpulse), Vdd(jpulse)
      real*8 Rxx(jpulse), Ryy(jpulse), phr(jplot), psv(jpulse)
     
      integer n, j, iprint, jp, istop, ipulse, ifill
      integer icycle_pulse, irep, nrep
      integer icycle_plot, idirn

      logical lstop, lpr

      intrinsic abs, atan, cos, sin, sqrt, acos

      external RK, indata_pha


c     constants
c     *********
      hpi = 2.0d00*atan(1.0d00)
      pi  = 2.0d00*hpi
      tpi = 4.0d00*hpi
      recip_root2=1.0d00/sqrt(2.0d00)

      open(unit=46,file='outdata_pha.txt',status='replace')
      open(unit=48,file='iterr_pha.txt',status='replace')
      write(48,984)
984   format('irep',9x,'time',2x,'Amplitude')

      istop=0
      call indata_pha(istop)
      ipulse = pulse_time_klys*fd
      ifill  = flip_time_klys*fd
      if(ipulse.gt.jpulse)then
        print*,'ipulse=',ipulse,' jpulse=',jpulse
        print*,'Increase array size jpulse'
        istop=1
      end if
      if(istop.eq.1) go to 3001

c     Set fill phase, flip phase and final phase from indata_pha.txt
c     **************************************************************
      psi0 = 0.0d00
      psi1 = phase_flip*pi
      if(psi1.gt.0.0d00)then
        psi2  = pi
        idirn = 1
      else
        psi2  =-pi
        idirn =-1
      end if

c     angular frequencies
c     *******************
      wd=tpi*fd   
      wcA=tpi*fcA
      wcB=tpi*fcB


c     solve envelop equations on every RF cycle
c     *****************************************
      dt     =  1.0d00/fd

c     assume flip is a time = 0 and filling starts at a negative time
c     **************************************************************
      tstart = -ifill*dt

      n=2*ifill

      open(file='results_pha.txt',unit=40,status='replace')

      write(40,930) 
930   format(10x,'time',6x,'cycle',3x,'    F_amp',3x,'  F_phase',
     +              3x,'    R_amp',3x,'  R_phase',3x,'      RAx',
     +              3x,'      RAy',3x,'      RBx',3x,'      RBy',
     +              3x,'      drr',3x,'      dri')

      hdt=dt*0.5d00
      dt6=dt/6.0d00
     
c     vary phase in small steps dpx and volts in steps dV
c     ***************************************************
c     dp1 is phase change per cycle at sudden phase switch
c     dp2 is phase change per cycle during ramp
c     dpp2 is to give a nonlinear ramp for phase

c     jp is determined by the klystron bandwidth
      fbandwidth = 20.0d00/12000.0d00
      dp1 = tpi*fbandwidth
      jp  = abs((psi1-psi0)/dp1)
      dp1 = (psi1-psi0)/jp
      dp2 = 0.0004d00*(psi2-psi1)

c     dV is determined by the klystron bandwidth
      dV = fbandwidth*Vd0
      Vd3= Vd0

c     loop to set initial phase ramp
c     ******************************
      do 3 j=1,jpulse
       if(j.lt.jp)then
        dpe(j)=dp1
       else
        dpe(j)=dp2
       end if
3     continue

      nrep = 200

      lstop=.false.


c     loop to correct phase ramp for flat top
c     ***************************************
      do 4 irep=1,nrep

      lpr=.true.

      t0=tstart
      ts=t0

      RAx=0.0d00
      RAy=0.0d00
      RBx=0.0d00
      RBy=0.0d00
      icycle_pulse = 0
      icycle_plot  = 0

      Vx = Vd0*cos(psi0)
      Vy = Vd0*sin(psi0)
      Vx = 0.0d00
      Vy = 0.0d00

      iprint = 5

c     loop to solve envelope equation with given phase ramp
c     *****************************************************
      do 1 j=1,n
        Vx0=Vx
        Vy0=Vy

        if(t0.ge.-20.0d-09) icycle_plot=icycle_plot+1

        if(t0.gt.0.0d00)then
          icycle_pulse = icycle_pulse + 1

          if(icycle_pulse.le.ipulse)then
            psi = psi + dpe(icycle_pulse)
            ph(icycle_pulse)=psi
          else
            psi = psi2
          end if

          if(icycle_pulse.le.ipulse)then
c           ****** in pulse voltage ******
            Vd = Vd0
            Vd3= Vd0
          else
c           ****** after pulse voltage run down and off *******
            if(Vd3.gt.0.0d00)then
              Vd3 = Vd3 - dV
              Vd  = Vd3
            else
              Vd = 0.0d00
            end if
          end if
        else
c         ****** filling voltage and phase ********
          icycle_pulse = 0
          psi=psi0
          Vd = Vd0
        end if
        psi_deg= psi*180.0d00/pi


c       set I and Q drive voltages and there differntials
c       *************************************************
        Vx = Vd*cos(psi)
        Vy = Vd*sin(psi)
        DVx=(Vx-Vx0)*fd
        DVy=(Vy-Vy0)*fd

c       step envelope equation for compressor cavity A
c       **********************************************
        g1=0.5*wd/QLA
        g2=wd-wcA
        g3=wcA/QEA-g1
        VAx = recip_root2*Vx
        VAy = recip_root2*Vy
        DVAx= recip_root2*DVx
        DVAy= recip_root2*DVy
        drr=-DVAx+g3*VAx-g2*VAy
        dri=-DVAy+g3*VAy+g2*VAx
        call RK(t0,hdt,dt6,RAx,RAy,g1,g2,drr,dri,t2)

c       step envelope equation for compressor cavity B
c       after applying 90 degree shift in hybrid coupler
c       ************************************************
        g1=0.5*wd/QLB
        g2=wd-wcB
        g3=wcB/QEB-g1
        VBx =-recip_root2*Vy
        VBy = recip_root2*Vx
        DVBx=-recip_root2*DVy
        DVBy= recip_root2*DVx
        drr=-DVBx+g3*VBx-g2*VBy
        dri=-DVBy+g3*VBy+g2*VBx
        call RK(t0,hdt,dt6,RBx,RBy,g1,g2,drr,dri,t2)


c       combine cavity outputs in Hybrid coupler
c       ****************************************
        Rx = recip_root2*(RBx - RAy)
        Ry = recip_root2*(RAx + RBy)
c       Rx+jRy is combined output after hybrid

        time = t2*1.0d09
      
        Ramp=sqrt(Rx**2+Ry**2)

        if(Ramp.gt.0.0d00)then
         if(Rx.ge.0.0d00)then
          Rpsi=asin(Ry/Ramp)
         else
          if(Ry.gt.0.0d00)then
           Rpsi= pi - asin(Ry/Ramp)
          else
           Rpsi=-pi - asin(Ry/Ramp)
          end if
         end if
        else
         Rpsi=0.0d00
        end if
        Rpsi_deg= Rpsi*180.0d00/pi


c       save output during compressed pulse 
c       ***********************************
        if((icycle_plot.gt.0).and.(icycle_plot.le.ipulse+1000))then
         tim(icycle_plot)  = time
         amp(icycle_plot)  = Ramp
         phr(icycle_plot)  = Rpsi_deg
        end if

        if((icycle_pulse.gt.0).and.(icycle_pulse.le.ipulse))then
         Rxx(icycle_pulse)  = Rx
         Ryy(icycle_pulse)  = Ry
         ampx(icycle_pulse) = Ramp
         psv(icycle_pulse)  = psi_deg
         Vdd(icycle_pulse)  = Vd
        end if

c       save every 5th output for final repetition
c       *******************************************
        if(iprint.eq.5)then
         if(lstop.or.(irep.eq.nrep))then
          v1=jp*drr/fd
          v2=jp*dri/fd
          write(40,902) time, icycle_pulse,Vd, psi_deg, Ramp, Rpsi_deg,
     +                  RAx, RAy, RBx, RBy, v1, v2
902       format(3x,f11.3,3x,i8,8(3x,f9.5),2(3x,1pe10.3))
         end if
         iprint=0
        end if
        iprint=iprint + 1

        if((t2.ge.errtime).and.lpr)then
          write(48,987)irep, time, Ramp
987       format(i4,2x,f11.3,2x,f9.5)
          lpr=.false.
        end if

        t0=t2
      
1     continue


c     adjust phase ramp for next repetion
c     ***********************************
      sum=0.0d00
c     fac1=0.021d00
c     fac2=0.0000025d00


      do 5 j=jp+1,ipulse-1

       if(j.eq.jp+1)then
         diff=ampx(j)-ampx(j+1)
       else
         diff=0.5d00*(ampx(j-1)-ampx(j+1))
       end if
       err=ampx(j)-ampx(jp)
     
       dpe(j)=dpe(j)+idirn*fac1*diff-idirn*fac2*err
       sum = sum + abs(diff)
5     continue
      dpe(ipulse)=2.0d00*dpe(ipulse-1)-1.0d00*dpe(ipulse-2)
      error=sum/(2.0d00*ipulse)

      if(irep.gt.1) error1=error
      write(*,966) error, irep
      write(46,966) error, irep
966   format('error=',1pe15.8,'  irep=',i3)


c     save selected inpulse output
c     ***************************
      do 6 j=1,ipulse + 1000
       if(irep.eq.1)  amp1(j)=amp(j)
       if(irep.eq.4)  amp2(j)=amp(j)
       if(irep.eq.8)  amp3(j)=amp(j)
       if(irep.eq.12)  amp4(j)=amp(j)
       if(irep.eq.1)  ph1(j)=phr(j)
       if(irep.eq.4)  ph2(j)=phr(j)
       if(irep.eq.8)  ph3(j)=phr(j)
       if(irep.eq.12)  ph4(j)=phr(j)
6     continue

      if(lstop)go to 7

      if(irep.gt.4)then
       if(error.lt.1.0e-8)then
        lstop=.true.
       end if
      end if

4     continue

7     close(unit=40,status='keep')

      open(file='pulse_pha.txt',unit=41,status='replace')
      write(41,988)
988   format(6x,'time',8x,'amp1',8x,'amp4',8x,'amp8',7x,'amp12',
     +       7x,'amp30',9x,'ph1',9x,'ph4',9x,'ph8',8x,'ph12',
     +       7x,'ph30',13X,'dp')
      do 2 j=1,ipulse + 1000
        write(41,977) tim(j),amp1(j),amp2(j),amp3(j),amp4(j),
     +                amp(j),ph1(j),ph2(j),ph3(j),ph4(j),
     +                phr(j),dpe(j)
977     format(11(f10.4,2x),2x,1pe13.5)
2     continue
      close(unit=41,status='keep')

      power_top = amp(1000)**2

      print*,'ipulse=',ipulse,'  power=',power_top
      write(46,*) 'ipulse=',ipulse,'  power=',power_top

      close(unit=48,status='keep')
3001  close(unit=46,status='keep')
       
      stop

      end


      subroutine RK(t0,hdt,dt6,AR0,AI0,f1,f2,outr,outi,t2)
c     ****************************************************
c     4th order Runge Kutta integration routine
      real*8 AR0, AR1, AR2, AR3, AR4, AI0, AI1, AI2, AI3, AI4
      real*8 DAR1, DAR2, DAR3, DAR4, DAI1, DAI2, DAI3, DAI4
      real*8 f1, f2, t0, t1, t2, dt6, hdt, outr, outi

      t1  = t0+hdt
      DAR1=-f1*AR0-f2*AI0+outr
      AR1 = AR0+hdt*DAR1
      DAI1=-f1*AI0+f2*AR0+outi
      AI1 = AI0+hdt*DAI1

      DAR2=-f1*AR1-f2*AI1+outr
      AR2 = AR0+hdt*DAR2
      DAI2=-f1*AI1+f2*AR1+outi
      AI2 = AI0+hdt*DAI2

      DAR3=-f1*AR2-f2*AI2+outr
      AR3 = AR0+2.0d00*hdt*DAR3
      DAI3=-f1*AI2+f2*AR2+outi
      AI3 = AI0+2.0d00*hdt*DAI3

      t2 = t1+hdt
      DAR4=-f1*AR3-f2*AI3+outr
      AR4 = AR0+dt6*(DAR1+2.0d00*(DAR2+DAR3)+DAR4)
      DAI4=-f1*AI3+f2*AR3+outi
      AI4 = AI0+dt6*(DAI1+2.0d00*(DAI2+DAI3)+DAI4)

      AR0=AR4
      AI0=AI4
      
      return
      end

      subroutine indata_pha(istop)
c     ****************************

      real*8 QLA, QEA, QLB, QEB, Q0A, Q0B
      real*8 fcA, fcB, fd, Vd0, fbandwidth
      real*8 pulse_time_mod, flip_time_klys, pulse_time_klys
      real*8 phase_flip, errtime, fac1, fac2

      common/com_imp/QLA, QEA, QLB, QEB, Q0A, Q0B,fcA, fcB, fd, Vd0,
     +               fbandwidth, pulse_time_mod, flip_time_klys,
     +               pulse_time_klys, phase_flip,
     +               errtime, fac1, fac2 


      integer*4 istop, nerror

      character*80 anything

      istop=0


      open(file='indata_pha.txt',unit=47,status='old') 
900   format(a)

c     cavity loaded and external Q factors
c     ************************************
c     Q0A=180000.0d00
c     Q0B=180000.0d00
c     QLA=40000.0d00
c     QLB=40000.0d00

c     drive frequency = fd, cavity A freq = fcA, cavity B freq = fcB
c     **************************************************************
c     fd  = 12.0d09
c     fcA = 12.0d09
c     fcB = 12.0d09

c     drive voltage
c     *************
c     Vd0 = 100.0d00
c     phase_flip=0.538d00

c     set fill time  and pulse length
c     *******************************
c     pulse_time_mod = 1.50d-6
c     flip_time_klys = 1.25d-6


      read(47,900) anything
      nerror=9521
      read(anything(41:),*,err=3000) Q0A
      write(*,9521) Q0A
      write(46,9521) Q0A
9521  format('Intrinsic Q factor cavity A             =',f12.1) 

      read(47,900) anything
      nerror=9522
      read(anything(41:),*,err=3000) Q0B
      write(*,9522) Q0B
      write(46,9522) Q0B
9522  format('Intrinsic Q factor cavity B             =',f12.1) 

      read(47,900) anything
      nerror=9523
      read(anything(41:),*,err=3000) QLA
      write(*,9523) QLA
      write(46,9523) QLA
9523  format('Loaded Q factor cavity A                =',f12.1) 

      read(47,900) anything
      nerror=9524
      read(anything(41:),*,err=3000) QLB
      write(*,9524) QLB
      write(46,9524) QLB
9524  format('Loaded Q factor cavity B                =',f12.1) 

      QEA=Q0A*QLA/(Q0A-QLA)
      QEB=Q0B*QLB/(Q0B-QLB)

      read(47,900) anything
      nerror=9525
      read(anything(41:),*,err=3000) fd
      write(*,9525) fd
      write(46,9525) fd
9525  format('Drive frequency                         =',f12.0) 

      read(47,900) anything
      nerror=9526
      read(anything(41:),*,err=3000) fcA
      write(*,9526) fcA
      write(46,9526) fcA
9526  format('Natural frequency of cavity A           =',f12.0) 

      read(47,900) anything
      nerror=9527
      read(anything(41:),*,err=3000) fcB
      write(*,9527) fcB
      write(46,9527) fcB
9527  format('Natural frequency of cavity B           =',f12.0) 

      read(47,900) anything
      nerror=9528
      read(anything(41:),*,err=3000) fbandwidth
      write(*,9528) fbandwidth
      write(46,9528) fbandwidth
9528  format('Klystron bandwidth                      =',f12.3) 

      read(47,900) anything
      nerror=9529
      read(anything(41:),*,err=3000) Vd0
      write(*,9529) Vd0
      write(46,9529) Vd0
9529  format('Drive voltage                           =',f12.3) 

      read(47,900) anything
      nerror=9530
      read(anything(41:),*,err=3000) phase_flip
      write(*,9530) phase_flip
      write(46,9530) phase_flip
9530  format('Phase change fraction at flip time      =',f12.6) 

      read(47,900) anything
      nerror=9531
      read(anything(41:),*,err=3000) pulse_time_mod
      write(*,9531) pulse_time_mod
      write(46,9531) pulse_time_mod
9531  format('Duration time of modulator pulse        =',1pe12.3) 

      read(47,900) anything
      nerror=9532
      read(anything(41:),*,err=3000) flip_time_klys
      write(*,9532) flip_time_klys
      write(46,9532) flip_time_klys
9532  format('Time for phase flip                     =',1pe12.3) 

      read(47,900) anything
      nerror=9531
      read(anything(41:),*,err=3000) fac1
      write(*,9533) fac1
      write(46,9533) fac1
9533  format('Proportional coefficent                 =',1pe12.3) 

      read(47,900) anything
      nerror=9532
      read(anything(41:),*,err=3000) fac2
      write(*,9532) fac2
      write(46,9532) fac2
9534  format('Integral coefficent                     =',1pe12.3) 

      pulse_time_klys= pulse_time_mod - flip_time_klys
      errtime        = 0.95d00*pulse_time_klys


      close(unit=47,status='keep')

      return

3000  istop=1

      print*,'nerror=',nerror
      close(unit=47,status='keep')
      return
      end