c...+....1....+....2....+....3....+....4....+....5....+....6....+....7....+
c     File:        mopagras.f (Fortran 77 source for the mopagras program)
c     Author:      Fredrik Jonsson <fj@phys.soton.ac.uk>
c     Date:        December 26, 2005
c     Last change: January 3, 2006
c     Description: The mopagras program calculates field envelopes of
c                  forward and backward travelling wave components inside
c                  a magneto-optically parametric bragg grating of sinusoi-
c                  dal spatial distribution of the refractive index.
c                     The program essentially wraps the calculating core
c                  provided by the mopagrat routine (see file mopagrat.f)
c                  into a stand-alone executable program, in which command-
c                  line parsing for parameters and output files is support-
c                  ed, as well as including the possibility to specifying
c                  different scan parameters, i.e. specifying different
c                  variables to be varied in investigating for example
c                  optimization or phase matching properties of the
c                  parametric process.
c
C     Copyright (C) 2006, Fredrik Jonsson <fj@phys.soton.ac.uk>
c...+....1....+....2....+....3....+....4....+....5....+....6....+....7....+
c     program main
      implicit logical (a-z)
      external zgaussj,mopagrat,ccubsolv
      logical verbose
      integer*4 j1,j2,k,numargs,oid_lcp,oid_rcp,num_scan_params
      integer*4 i_alpha,num_alpha_steps
      integer*4 i_beta,num_beta_steps
      integer*4 i_gamma,num_gamma_steps
      integer*4 i_delta,num_delta_steps
      integer*4 i_eta,num_eta_steps
      integer*4 i_kappa,num_kappa_steps
      integer*4 i_length,num_length_steps
      integer*4 i_zeta,num_zeta_steps
      integer*4 i_magfield,num_magfield_steps
      real*8 alpha,alphastart,alphastop,dalpha
      real*8 beta,betastart,betastop,dbeta
      real*8 gamma,gammastart,gammastop,dgamma
      real*8 delta,deltastart,deltastop,ddelta
      real*8 eta,etastart,etastop,deta
      real*8 kappa,kappastart,kappastop,dkappa
      real*8 length,lengthstart,lengthstop,dlength
      real*8 zeta,zetastart,zetastop,dzeta
      real*8 magfield,magfieldstart,magfieldstop,dmagfield
      real*8 tmpreal,tmpgamma,tmpdelta
      complex*8 a1iplus,a2iplus,a1fplus,a1bplus,a2fplus
      complex*8 a1iminus,a2iminus,a1fminus,a1bminus,a2fminus
      real*8 a1famplplus,a1bamplplus,a2famplplus
      real*8 a1famplminus,a1bamplminus,a2famplminus
      character*48 arg,outfilename_lcp,outfilename_rcp

c     *********************************************************************
c     Initialization of variables with default values
c     *********************************************************************
      verbose=false
      outfilename_lcp='out.lcp.dat'
      outfilename_rcp='out.rcp.dat'
      oid_lcp=10
      oid_rcp=11
      num_scan_params=0
      a1fplus=cmplx(0.0,0.0)
      a1bplus=cmplx(0.0,0.0)
      a2fplus=cmplx(0.0,0.0)
      a1fminus=cmplx(0.0,0.0)
      a1bminus=cmplx(0.0,0.0)
      a2fminus=cmplx(0.0,0.0)
      alpha=0.0d0
      alphastart=0.0d0
      alphastop=0.0d0
      dalpha=-1.0d0
      num_alpha_steps=1
      beta=0.0d0
      betastart=0.0d0
      betastop=1.0d0
      dbeta=-1.0d0
      num_beta_steps=1
      gamma=0.0d0
      gammastart=0.0d0
      gammastop=1.0d0
      dgamma=-1.0d0
      num_gamma_steps=1
      delta=0.0d0
      deltastart=0.0d0
      deltastop=1.0d0
      ddelta=-1.0d0
      num_delta_steps=1
      eta=0.0d0
      etastart=0.0d0
      etastop=1.0d0
      deta=-1.0d0
      num_eta_steps=1
      kappa=0.0d0
      kappastart=0.0d0
      kappastop=1.0d0
      dkappa=-1.0d0
      num_kappa_steps=1
      length=0.0d0
      lengthstart=0.0d0
      lengthstop=1.0d0
      dlength=-1.0d0
      num_length_steps=1
      zeta=0.0d0
      zetastart=0.0d0
      zetastop=1.0d0
      dzeta=-1.0d0
      num_zeta_steps=1
      magfield=0.0d0
      magfieldstart=0.0d0
      magfieldstop=0.0d0
      dmagfield=1.0d0
      num_magfield_steps=1
      a1iplus=cmplx(0.0,0.0)
      a2iplus=cmplx(0.0,0.0)
      a1iminus=cmplx(0.0,0.0)
      a2iminus=cmplx(0.0,0.0)
      a1famplplus=cmplx(0.0,0.0)
      a1bamplplus=cmplx(0.0,0.0)
      a2famplplus=cmplx(0.0,0.0)
      a1famplminus=cmplx(0.0,0.0)
      a1bamplminus=cmplx(0.0,0.0)
      a2famplminus=cmplx(0.0,0.0)

c     *********************************************************************
c     Get the input and output file names from the argument list.
c     in parsing the list of numerical intervals of the form
c
c       <option switch> <start value>:<stop value>:<# scan points>
c
c     The integers j1 and j2 are here used to keep track of the positions
c     of the colons in the parsed arg string, providing the means of a
c     compact and consistent notation in the command line syntax.
c     *********************************************************************
      numargs = iargc()
      k=1
      do while(k.le.numargs)
         call getarg(k,arg)
         k=k+1
         if (arg.eq.'--verbose') then
            verbose=true
            write(6,10) 'entering verbose mode'
 10         format('mopagras: ',a)
         elseif ((arg.eq.'-olcp').or.(arg.eq.'--outputfile_lcp')) then
            call getarg(k,outfilename_lcp)
            k=k+1
         elseif ((arg.eq.'-orcp').or.(arg.eq.'--outputfile_rcp')) then
            call getarg(k,outfilename_rcp)
            k=k+1
         elseif ((arg.eq.'--a1iplus').or.(arg.eq.'--signal_lcp')) then
            call getarg(k,arg)
            k=k+1
            read(arg,*) tmpreal
            a1iplus=cmplx(tmpreal,0.0)
         elseif ((arg.eq.'--a2iplus').or.(arg.eq.'--idler_lcp')) then
            call getarg(k,arg)
            k=k+1
            read(arg,*) tmpreal
            a2iplus=cmplx(tmpreal,0.0)
         elseif ((arg.eq.'--a1iminus').or.(arg.eq.'--signal_rcp')) then
            call getarg(k,arg)
            k=k+1
            read(arg,*) tmpreal
            a1iminus=cmplx(tmpreal,0.0)
         elseif ((arg.eq.'--a2iminus').or.(arg.eq.'--idler_rcp')) then
            call getarg(k,arg)
            k=k+1
            read(arg,*) tmpreal
            a2iminus=cmplx(tmpreal,0.0)
         elseif (arg.eq.'--magfield') then
            call getarg(k,arg)
            k=k+1
            j1=index(arg,':')
            read(arg(1:j1-1),*) magfieldstart
            arg(1:j1)=' '
            j2=index(arg,':')
            read(arg(j1+1:j2-1),*) magfieldstop
            read(arg(j2+1:),*) num_magfield_steps
         elseif (arg.eq.'--alpha') then
            call getarg(k,arg)
            k=k+1
            j1=index(arg,':')
            read(arg(1:j1-1),*) alphastart
            arg(1:j1)=' '
            j2=index(arg,':')
            read(arg(j1+1:j2-1),*) alphastop
            read(arg(j2+1:),*) num_alpha_steps
         elseif (arg.eq.'--beta') then
            call getarg(k,arg)
            k=k+1
            j1=index(arg,':')
            read(arg(1:j1-1),*) betastart
            arg(1:j1)=' '
            j2=index(arg,':')
            read(arg(j1+1:j2-1),*) betastop
            read(arg(j2+1:),*) num_beta_steps
         elseif (arg.eq.'--gamma') then
            call getarg(k,arg)
            k=k+1
            j1=index(arg,':')
            read(arg(1:j1-1),*) gammastart
            arg(1:j1)=' '
            j2=index(arg,':')
            read(arg(j1+1:j2-1),*) gammastop
            read(arg(j2+1:),*) num_gamma_steps
         elseif (arg.eq.'--delta') then
            call getarg(k,arg)
            k=k+1
            j1=index(arg,':')
            read(arg(1:j1-1),*) deltastart
            arg(1:j1)=' '
            j2=index(arg,':')
            read(arg(j1+1:j2-1),*) deltastop
            read(arg(j2+1:),*) num_delta_steps
         elseif (arg.eq.'--eta') then
            call getarg(k,arg)
            k=k+1
            j1=index(arg,':')
            read(arg(1:j1-1),*) etastart
            arg(1:j1)=' '
            j2=index(arg,':')
            read(arg(j1+1:j2-1),*) etastop
            read(arg(j2+1:),*) num_eta_steps
         elseif (arg.eq.'--kappa') then
            call getarg(k,arg)
            k=k+1
            j1=index(arg,':')
            read(arg(1:j1-1),*) kappastart
            arg(1:j1)=' '
            j2=index(arg,':')
            read(arg(j1+1:j2-1),*) kappastop
            read(arg(j2+1:),*) num_kappa_steps
         elseif (arg.eq.'--length') then
            call getarg(k,arg)
            k=k+1
            j1=index(arg,':')
            read(arg(1:j1-1),*) lengthstart
            arg(1:j1)=' '
            j2=index(arg,':')
            read(arg(j1+1:j2-1),*) lengthstop
            read(arg(j2+1:),*) num_length_steps
         elseif (arg.eq.'--zeta') then
            call getarg(k,arg)
            k=k+1
            j1=index(arg,':')
            read(arg(1:j1-1),*) zetastart
            arg(1:j1)=' '
            j2=index(arg,':')
            read(arg(j1+1:j2-1),*) zetastop
            read(arg(j2+1:),*) num_zeta_steps
         else
            write(6,20) 'Error: Unrecognized command line option: ',arg
 20         format('mopagras: ',a41,a26)
            stop
         endif
      enddo

c     *********************************************************************
c     Initialize the step lengths to be taken for each scan parameter.
c     *********************************************************************
      if (num_magfield_steps.gt.1) then
         dmagfield=(magfieldstop-magfieldstart)/(num_magfield_steps-1)
         num_scan_params=num_scan_params+1
      else
         dzeta=(zetastop-zetastart)
      endif
      if (num_alpha_steps.gt.1) then
         dalpha=(alphastop-alphastart)/(num_alpha_steps-1)
         num_scan_params=num_scan_params+1
      else
         dalpha=(alphastop-alphastart)
      endif
      if (num_beta_steps.gt.1) then
         dbeta=(betastop-betastart)/(num_beta_steps-1)
         num_scan_params=num_scan_params+1
      else
         dbeta=(betastop-betastart)
      endif
      if (num_gamma_steps.gt.1) then
         dgamma=(gammastop-gammastart)/(num_gamma_steps-1)
         num_scan_params=num_scan_params+1
      else
         dgamma=(gammastop-gammastart)
      endif
      if (num_delta_steps.gt.1) then
         ddelta=(deltastop-deltastart)/(num_delta_steps-1)
         num_scan_params=num_scan_params+1
      else
         ddelta=(deltastop-deltastart)
      endif
      if (num_eta_steps.gt.1) then
         deta=(etastop-etastart)/(num_eta_steps-1)
         num_scan_params=num_scan_params+1
      else
         deta=(etastop-etastart)
      endif
      if (num_kappa_steps.gt.1) then
         dkappa=(kappastop-kappastart)/(num_kappa_steps-1)
         num_scan_params=num_scan_params+1
      else
         dkappa=(kappastop-kappastart)
      endif
      if (num_length_steps.gt.1) then
         dlength=(lengthstop-lengthstart)/(num_length_steps-1)
         num_scan_params=num_scan_params+1
      else
         dlength=(lengthstop-lengthstart)
      endif
      if (num_zeta_steps.gt.1) then
         dzeta=(zetastop-zetastart)/(num_zeta_steps-1)
         num_scan_params=num_scan_params+1
      else
         dzeta=(zetastop-zetastart)
      endif
      if (num_scan_params.gt.1) then
         write(6,20) 'error: more than one scan parameter specified!'
         write(6,20) '(please check parameter intervals for errors)'
         stop
      endif

c     *********************************************************************
c     If in verbose mode, show the parsed scan parameter interals at
c     standard terminal output.
c     *********************************************************************
      if (verbose.eqv.true) then
         write(6,30) '------------------------------------------------'
 30      format('------------',a,'------------')
         write(6,40) 'Summary of parameters parsed from command line'
 40      format(a)
         write(6,50) numargs
 50      format('number of arguments parsed: ',i2)
         write(6,60) 'Input signal field amplitude (lcp): ',a1iplus
         write(6,60) 'Input signal field amplitude (rcp): ',a1iminus
         write(6,60) 'Input idler field amplitude (lcp):  ',a2iplus
         write(6,60) 'Input idler field amplitude (rcp):  ',a2iminus
 60      format(3x,a36,e15.7,e15.7)
         write(6,70) 'Start value','Stop value','Increment','Scanning'
 70      format(10x,a16,a16,a16,a13)
         arg='(no scan)'
         if (num_magfield_steps.gt.1) arg='(scan)'
         write(6,80) 'magfield: ',magfieldstart,magfieldstop,
     +        dmagfield,arg
         arg='(no scan)'
         if (num_alpha_steps.gt.1) arg='(scan)'
         write(6,80) 'alpha:    ',alphastart,alphastop,dalpha,arg
         arg='(no scan)'
         if (num_beta_steps.gt.1) arg='(scan)'
         write(6,80) 'beta:     ',betastart,betastop,dbeta,arg
         arg='(no scan)'
         if (num_gamma_steps.gt.1) arg='(scan)'
         write(6,80) 'gamma:    ',gammastart,gammastop,dgamma,arg
         arg='(no scan)'
         if (num_delta_steps.gt.1) arg='(scan)'
         write(6,80) 'delta:    ',deltastart,deltastop,ddelta,arg
         arg='(no scan)'
         if (num_eta_steps.gt.1) arg='(scan)'
         write(6,80) 'eta:      ',etastart,etastop,deta,arg
         arg='(no scan)'
         if (num_kappa_steps.gt.1) arg='(scan)'
         write(6,80) 'kappa:    ',kappastart,kappastop,dkappa,arg
         arg='(no scan)'
         if (num_length_steps.gt.1) arg='(scan)'
         write(6,80) 'length:   ',lengthstart,lengthstop,dlength,arg
         arg='(no scan)'
         if (num_zeta_steps.gt.1) arg='(scan)'
         write(6,80) 'zeta:     ',zetastart,zetastop,dzeta,arg
 80      format(3x,a10,e15.7,'  ',e15.7,'  ',e15.7,' ',a9)
         write(6,30) '------------------------------------------------'
      endif

c     *********************************************************************
c     Open the file specified by the string ouputfilename for output.
c     *********************************************************************
      open(oid_lcp,file=outfilename_lcp)
      open(oid_rcp,file=outfilename_rcp)

c     *********************************************************************
c     Enter the main loop of parameter scanning in the calculation of
c     grating field envelopes.
c     *********************************************************************
      do 100,i_alpha=1,num_alpha_steps
         alpha=alphastart+(i_alpha-1)*dalpha
         do 110,i_beta=1,num_beta_steps
            beta=betastart+(i_beta-1)*dbeta
            do 120,i_gamma=1,num_gamma_steps
               gamma=gammastart+(i_gamma-1)*dgamma
               do 130,i_delta=1,num_delta_steps
                  delta=deltastart+(i_delta-1)*ddelta
                  do 140,i_eta=1,num_eta_steps
                     eta=etastart+(i_eta-1)*deta
                     do 150,i_kappa=1,num_kappa_steps
                        kappa=kappastart+(i_kappa-1)*dkappa
                        do 160,i_length=1,num_length_steps
                           length=lengthstart+(i_length-1)*dlength
                           do 170,i_zeta=1,num_zeta_steps
                              zeta=zetastart+(i_zeta-1)*dzeta
                              do 180,i_magfield=1,num_magfield_steps
                                 magfield=magfieldstart
     +                                +(i_magfield-1)*dmagfield
                           if (num_magfield_steps.gt.1) then
                              tmpgamma=magfield*gamma
                              tmpdelta=magfield*delta
                           call mopagrat(alpha,beta,tmpgamma,
     +                       tmpdelta,eta,kappa,length,zeta,
     +                       a1iplus,a2iplus,a1iminus,a2iminus,
     +                       a1famplplus,a1bamplplus,a2famplplus,
     +                       a1famplminus,a1bamplminus,a2famplminus,
     +                       a1fplus,a1bplus,a2fplus,
     +                       a1fminus,a1bminus,a2fminus)
                           else
                           call mopagrat(alpha,beta,gamma,
     +                       delta,eta,kappa,length,zeta,
     +                       a1iplus,a2iplus,a1iminus,a2iminus,
     +                       a1famplplus,a1bamplplus,a2famplplus,
     +                       a1famplminus,a1bamplminus,a2famplminus,
     +                       a1fplus,a1bplus,a2fplus,
     +                       a1fminus,a1bminus,a2fminus)
                           endif

                           if (num_alpha_steps.gt.1) then
                              write(oid_lcp,90) alpha,
     +                        a1famplplus,a1bamplplus,a2famplplus
                              write(oid_rcp,90) alpha,
     +                        a1famplminus,a1bamplminus,a2famplminus
                           endif
                           if (num_beta_steps.gt.1) then
                              write(oid_lcp,90) beta,
     +                        a1famplplus,a1bamplplus,a2famplplus
                              write(oid_rcp,90) beta,
     +                        a1famplminus,a1bamplminus,a2famplminus
                           endif
                           if (num_magfield_steps.gt.1) then
                              write(oid_lcp,90) magfield,
     +                        a1famplplus,a1bamplplus,a2famplplus
                              write(oid_rcp,90) magfield,
     +                        a1famplminus,a1bamplminus,a2famplminus
                           else
                              if (num_gamma_steps.gt.1) then
                                 write(oid_lcp,90) gamma,
     +                        a1famplplus,a1bamplplus,a2famplplus
                                 write(oid_rcp,90) gamma,
     +                        a1famplminus,a1bamplminus,a2famplminus
                              endif
                              if (num_delta_steps.gt.1) then
                                 write(oid_lcp,90) delta,
     +                        a1famplplus,a1bamplplus,a2famplplus
                                 write(oid_rcp,90) delta,
     +                        a1famplminus,a1bamplminus,a2famplminus
                              endif
                           endif
                           if (num_eta_steps.gt.1) then
                              write(oid_lcp,90) eta,
     +                        a1famplplus,a1bamplplus,a2famplplus
                              write(oid_rcp,90) eta,
     +                        a1famplminus,a1bamplminus,a2famplminus
                           endif
                           if (num_kappa_steps.gt.1) then
                              write(oid_lcp,90) kappa,
     +                        a1famplplus,a1bamplplus,a2famplplus
                              write(oid_rcp,90) kappa,
     +                        a1famplminus,a1bamplminus,a2famplminus
                           endif
                           if (num_length_steps.gt.1) then
                              write(oid_lcp,90) length,
     +                        a1famplplus,a1bamplplus,a2famplplus
                              write(oid_rcp,90) length,
     +                        a1famplminus,a1bamplminus,a2famplminus
                           endif
                           if (num_zeta_steps.gt.1) then
                              write(oid_lcp,90) zeta,
     +                        a1famplplus,a1bamplplus,a2famplplus
                              write(oid_rcp,90) zeta,
     +                        a1famplminus,a1bamplminus,a2famplminus
                           endif
 180                          continue
 170                       continue
 160                    continue
 150                 continue
 140              continue
 130           continue
 120        continue
 110     continue
 100  continue
 90   format(1p,4(1x,e16.8))

c     *********************************************************************
c     Close the output files.
c     *********************************************************************
      close(oid_lcp)
      close(oid_rcp)

      stop
      end