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 has variable I and Q inputs
c     determined by a PI controller that seek constant amplitude and phase for the 
c     compressed pulse on the output.

c     To complile use ANSI FORTRAN77 
c     To run the datafile "indata_amp.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 Vd0x, Vd1x, Vd2x, Vd0y, Vd1y, Vd2y 
      real*8 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, Vd0x,Vd1x,Vd2x, Vd0y,Vd1y,Vd2y,
     +               errtime, fac1, fac2 

      real*8 hdt, dt6, dt
      real*8 hpi, pi, tpi, recip_root2
      real*8 Vd1, Vd2, Vd, v1, v2
      real*8 wcA, wcB, wd
      real*8 t0, t2, tstart, ts, time
      real*8 psi2, psi1, psi0, 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 DVD1, DVD2
      real*8 g3, Ramp, Rx, RY
      real*8 vmax
      real*8 err, diff, error, error1, sum
      real*8 power_top

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

      logical lstop, lpr

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

      external RK, indata_amp

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_amp.txt',status='replace')
      open(unit=48,file='iterr.txt',status='replace')
      write(48,984)
984   format('irep',9x,'time',2x,'Amplitude')

      istop=0
      call indata_amp(istop)
      ipulse = pulse_time_klys*fd
      ifill  = flip_time_klys*fd
      if((ipulse-500).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 IQ, flip IQ and final IQ from indata_amp.txt
c     *****************************************************
c     vd is the prefix for the klystron voltage 
c     start at (vd0x, vd0y)

c     move rapidly to (vd1x, vd1y) with unifrom steps in amplitude and phase
c     as permitted by bandwidth

c     move new point intially given by (vd2x, vd2y) so as to keep certain output
c     parameters derived from (Rx, Ry) such as amplitude or phase constant.

      vmax=sqrt(Vd0x**2+Vd0y**2)

      Vd0  = sqrt(Vd0x**2+Vd0y**2)
      Vd1  = sqrt(Vd1x**2+Vd1y**2)
      Vd2  = sqrt(Vd2x**2+Vd2y**2)

      psi0=0.0d00
      if(vd0.gt.1.0e-20) psi0 = acos(Vd0x/Vd0)

      psi1=0.0d00
      if(vd1.gt.1.0e-20) psi1 = acos(Vd1x/Vd1)

      psi2=0.0d00
      if(vd2.gt.1.0e-20) psi2 = acos(Vd2x/Vd2)

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_amp.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     dVd1 is amplitude 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     = (psi1-psi0)/dp1
      jph=0.85d00*jp
      if(jp.lt.4)  jp=4
      if(jph.lt.2) jph=3
      dp1    = (psi1-psi0)/(jp-1)
      dVd1   = (Vd1-Vd0)/(jph-1)

      dp2    = 0.0004d00*(psi2-psi1)
      dVd2   = 0.0004d00*(Vd2-Vd1)

c     loop to set initial phase ramp
c     ******************************
      Vx=Vd0x
      Vy=Vd0y
      
      Vd    = sqrt(Vx**2+Vy**2)
      if(Vd.gt.0.0d00)then
       if(Vx.ge.0.0d00)then
        psi=asin(Vy/Vd)
       else
        if(Vy.gt.0.0d00)then
         psi= pi - asin(Vy/Vd)
        else
         psi=-pi - asin(Vy/Vd)
        end if
       end if
      else
       psi=0.0d00
      end if
      psi_deg= psi*180.0d00/pi
      
      do 3 j=1,jpulse
       if(j.lt.jp)then
        dpe(j)=dp1
        if(j.lt.jph)then          
          dVe(j)=dVd1
        else
          dVe(j)=0.0d00
        end if
       else
        dpe(j)=dp2
        dVe(j)=dVd2
       end if
       psi = psi +dpe(j)
       Vd  = Vd + dVe(j)
       Vx0 = Vx
       Vy0 = Vy
       Vx = Vd*cos(psi)
       Vy = Vd*sin(psi)
       dVdx(j) = Vx - Vx0
       dVdy(j) = Vy - Vy0
3     continue

      nrep = 24

      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 = Vd0x
      Vy = Vd0y
      Vx0= 0.0d00
      Vy0= 0.0d00

      iprint = 5

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

        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
            Vx = Vx + dVdx(icycle_pulse)
            Vy = Vy + dVdy(icycle_pulse)
          else
            if(abs(Vx).gt.0.0d00)  Vx=0.95d00*Vx
            if(abs(Vy).gt.0.0d00)  Vy=0.95d00*Vy
          end if
        else
c         ****** filling voltage and phase ********
          icycle_pulse = 0
          Vx=Vd0x
          Vy=Vd0y
        end if

        Vd    = sqrt(Vx**2+Vy**2)
        if(Vd.gt.Vmax)then
          Vx=Vmax*cos(psi)
          Vy=Vmax*sin(psi)
          dVdx(icycle_pulse)=Vx-Vx0
          dVdy(icycle_pulse)=Vy-Vy0
        end if

        if(Vd.gt.0.0d00)then
         if(Vx.ge.0.0d00)then
          psi=asin(Vy/Vd)
         else
          if(Vy.gt.0.0d00)then
           psi= pi - asin(Vy/Vd)
          else
           psi=-pi - asin(Vy/Vd)
          end if
         end if
        else
         psi=0.0d00
        end if
        psi_deg= psi*180.0d00/pi


c       set I and Q drive voltages and there differentials
c       *************************************************
        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
         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

        Vx0=Vx
        Vy0=Vy
      
1     continue


c     adjust phase ramp for next repetion
c     ***********************************
      sum=0.0d00
c     fac1=0.5d00
c     fac2=0.0013d00

      do 5 j=jp-1,ipulse-1

       if(j.eq.jp-1)then
         diff=Rxx(j)-Rxx(j+1)
       else
         diff=0.5d00*(Rxx(j-1)-Rxx(j+1))
       end if
       err=Rxx(j)-Rxx(jp)
       dVdy(j)=dVdy(j) + fac1*diff-fac2*err
       sum = sum + abs(diff)

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

      if(irep.gt.1) error1=error
      print*,'err=',error,'  irep=',irep
  
c     save selected inpulse output
c     ***************************
      do 6 j=1,ipulse + 1000
       if(irep.eq.1)  amp1(j)=amp(j)
       if(irep.eq.2)  amp2(j)=amp(j)
       if(irep.eq.3)  amp3(j)=amp(j)
       if(irep.eq.4)  amp4(j)=amp(j)
       if(irep.eq.5)  amp5(j)=amp(j)      
       if(irep.eq.1)  ph1(j)=phr(j)
       if(irep.eq.2)  ph2(j)=phr(j)
       if(irep.eq.3)  ph3(j)=phr(j)
       if(irep.eq.4)  ph4(j)=phr(j)
       if(irep.eq.5)  ph5(j)=phr(j)
6     continue

      if(lstop)go to 7

      if(irep.gt.5)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_amp.txt',unit=41,status='replace')
      write(41,988)
988   format(6x,'time',8x,'amp1',8x,'amp2',8x,'amp3',8x,'amp4',
     +       8x,'amp5',7x,'amp30',9x,'ph1',9x,'ph2',9x,'ph3',
     +       9x,'ph4',9x,'ph5',8x,'ph30',13x,'dp')
      do 2 j=1,ipulse+1000
        write(41,977) tim(j),amp1(j),amp2(j),amp3(j),amp4(j),
     +                amp5(j),amp(j),ph1(j),ph2(j),ph3(j),
     +                ph4(j),ph5(j),phr(j),dpe(j)
977     format(13(f10.4,2x),1pe13.5,2x)
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     ****************************************************
      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_amp(istop)
c     ****************************
c23456789012345678901234567890123456789012345678901234567890123456789012
      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 Vd0x, Vd1x, Vd2x, Vd0y, Vd1y, Vd2y 
      real*8 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, Vd0x,Vd1x,Vd2x, Vd0y,Vd1y,Vd2y,
     +               errtime, fac1, fac2 


      integer*4 istop, nerror

      character*80 anything

      istop=0




      open(file='indata_amp.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


      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) Vd0x, Vd1x, Vd2x
      write(*,9529)  Vd0x, Vd1x, Vd2x
      write(46,9529)  Vd0x, Vd1x, Vd2x
9529  format('In phase Drive Voltages                 =',3(f12.3,2x)) 

      read(47,900) anything
      nerror=9530
      read(anything(41:),*,err=3000) Vd0y, Vd1y, Vd2y
      write(*,9530)  Vd0y, Vd1y, Vd2y
      write(46,9530)  Vd0y, Vd1y, Vd2y
9530  format('Quadrature Drive Voltages               =',3(f12.3,2x)) 

      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