Search:

Return to previous page

Contents of file 'mopagras.f':



    1   c...+....1....+....2....+....3....+....4....+....5....+....6....+....7....+
    2   c     File:        mopagras.f (Fortran 77 source for the mopagras program)
    3   c     Author:      Fredrik Jonsson <fj@phys.soton.ac.uk>
    4   c     Date:        December 26, 2005
    5   c     Last change: January 3, 2006
    6   c     Description: The mopagras program calculates field envelopes of
    7   c                  forward and backward travelling wave components inside
    8   c                  a magneto-optically parametric bragg grating of sinusoi-
    9   c                  dal spatial distribution of the refractive index.
   10   c                     The program essentially wraps the calculating core
   11   c                  provided by the mopagrat routine (see file mopagrat.f)
   12   c                  into a stand-alone executable program, in which command-
   13   c                  line parsing for parameters and output files is support-
   14   c                  ed, as well as including the possibility to specifying
   15   c                  different scan parameters, i.e. specifying different
   16   c                  variables to be varied in investigating for example
   17   c                  optimization or phase matching properties of the
   18   c                  parametric process.
   19   c
   20   C     Copyright (C) 2006, Fredrik Jonsson <fj@phys.soton.ac.uk>
   21   c...+....1....+....2....+....3....+....4....+....5....+....6....+....7....+
   22   c     program main
   23         implicit logical (a-z)
   24         external zgaussj,mopagrat,ccubsolv
   25         logical verbose
   26         integer*4 j1,j2,k,numargs,oid_lcp,oid_rcp,num_scan_params
   27         integer*4 i_alpha,num_alpha_steps
   28         integer*4 i_beta,num_beta_steps
   29         integer*4 i_gamma,num_gamma_steps
   30         integer*4 i_delta,num_delta_steps
   31         integer*4 i_eta,num_eta_steps
   32         integer*4 i_kappa,num_kappa_steps
   33         integer*4 i_length,num_length_steps
   34         integer*4 i_zeta,num_zeta_steps
   35         integer*4 i_magfield,num_magfield_steps
   36         real*8 alpha,alphastart,alphastop,dalpha
   37         real*8 beta,betastart,betastop,dbeta
   38         real*8 gamma,gammastart,gammastop,dgamma
   39         real*8 delta,deltastart,deltastop,ddelta
   40         real*8 eta,etastart,etastop,deta
   41         real*8 kappa,kappastart,kappastop,dkappa
   42         real*8 length,lengthstart,lengthstop,dlength
   43         real*8 zeta,zetastart,zetastop,dzeta
   44         real*8 magfield,magfieldstart,magfieldstop,dmagfield
   45         real*8 tmpreal,tmpgamma,tmpdelta
   46         complex*8 a1iplus,a2iplus,a1fplus,a1bplus,a2fplus
   47         complex*8 a1iminus,a2iminus,a1fminus,a1bminus,a2fminus
   48         real*8 a1famplplus,a1bamplplus,a2famplplus
   49         real*8 a1famplminus,a1bamplminus,a2famplminus
   50         character*48 arg,outfilename_lcp,outfilename_rcp
   51   
   52   c     *********************************************************************
   53   c     Initialization of variables with default values
   54   c     *********************************************************************
   55         verbose=false
   56         outfilename_lcp='out.lcp.dat'
   57         outfilename_rcp='out.rcp.dat'
   58         oid_lcp=10
   59         oid_rcp=11
   60         num_scan_params=0
   61         a1fplus=cmplx(0.0,0.0)
   62         a1bplus=cmplx(0.0,0.0)
   63         a2fplus=cmplx(0.0,0.0)
   64         a1fminus=cmplx(0.0,0.0)
   65         a1bminus=cmplx(0.0,0.0)
   66         a2fminus=cmplx(0.0,0.0)
   67         alpha=0.0d0
   68         alphastart=0.0d0
   69         alphastop=0.0d0
   70         dalpha=-1.0d0
   71         num_alpha_steps=1
   72         beta=0.0d0
   73         betastart=0.0d0
   74         betastop=1.0d0
   75         dbeta=-1.0d0
   76         num_beta_steps=1
   77         gamma=0.0d0
   78         gammastart=0.0d0
   79         gammastop=1.0d0
   80         dgamma=-1.0d0
   81         num_gamma_steps=1
   82         delta=0.0d0
   83         deltastart=0.0d0
   84         deltastop=1.0d0
   85         ddelta=-1.0d0
   86         num_delta_steps=1
   87         eta=0.0d0
   88         etastart=0.0d0
   89         etastop=1.0d0
   90         deta=-1.0d0
   91         num_eta_steps=1
   92         kappa=0.0d0
   93         kappastart=0.0d0
   94         kappastop=1.0d0
   95         dkappa=-1.0d0
   96         num_kappa_steps=1
   97         length=0.0d0
   98         lengthstart=0.0d0
   99         lengthstop=1.0d0
  100         dlength=-1.0d0
  101         num_length_steps=1
  102         zeta=0.0d0
  103         zetastart=0.0d0
  104         zetastop=1.0d0
  105         dzeta=-1.0d0
  106         num_zeta_steps=1
  107         magfield=0.0d0
  108         magfieldstart=0.0d0
  109         magfieldstop=0.0d0
  110         dmagfield=1.0d0
  111         num_magfield_steps=1
  112         a1iplus=cmplx(0.0,0.0)
  113         a2iplus=cmplx(0.0,0.0)
  114         a1iminus=cmplx(0.0,0.0)
  115         a2iminus=cmplx(0.0,0.0)
  116         a1famplplus=cmplx(0.0,0.0)
  117         a1bamplplus=cmplx(0.0,0.0)
  118         a2famplplus=cmplx(0.0,0.0)
  119         a1famplminus=cmplx(0.0,0.0)
  120         a1bamplminus=cmplx(0.0,0.0)
  121         a2famplminus=cmplx(0.0,0.0)
  122   
  123   c     *********************************************************************
  124   c     Get the input and output file names from the argument list.
  125   c     in parsing the list of numerical intervals of the form
  126   c
  127   c       <option switch> <start value>:<stop value>:<# scan points>
  128   c
  129   c     The integers j1 and j2 are here used to keep track of the positions
  130   c     of the colons in the parsed arg string, providing the means of a
  131   c     compact and consistent notation in the command line syntax.
  132   c     *********************************************************************
  133         numargs = iargc()
  134         k=1
  135         do while(k.le.numargs)
  136            call getarg(k,arg)
  137            k=k+1
  138            if (arg.eq.'--verbose') then
  139               verbose=true
  140               write(6,10) 'entering verbose mode'
  141    10         format('mopagras: ',a)
  142            elseif ((arg.eq.'-olcp').or.(arg.eq.'--outputfile_lcp')) then
  143               call getarg(k,outfilename_lcp)
  144               k=k+1
  145            elseif ((arg.eq.'-orcp').or.(arg.eq.'--outputfile_rcp')) then
  146               call getarg(k,outfilename_rcp)
  147               k=k+1
  148            elseif ((arg.eq.'--a1iplus').or.(arg.eq.'--signal_lcp')) then
  149               call getarg(k,arg)
  150               k=k+1
  151               read(arg,*) tmpreal
  152               a1iplus=cmplx(tmpreal,0.0)
  153            elseif ((arg.eq.'--a2iplus').or.(arg.eq.'--idler_lcp')) then
  154               call getarg(k,arg)
  155               k=k+1
  156               read(arg,*) tmpreal
  157               a2iplus=cmplx(tmpreal,0.0)
  158            elseif ((arg.eq.'--a1iminus').or.(arg.eq.'--signal_rcp')) then
  159               call getarg(k,arg)
  160               k=k+1
  161               read(arg,*) tmpreal
  162               a1iminus=cmplx(tmpreal,0.0)
  163            elseif ((arg.eq.'--a2iminus').or.(arg.eq.'--idler_rcp')) then
  164               call getarg(k,arg)
  165               k=k+1
  166               read(arg,*) tmpreal
  167               a2iminus=cmplx(tmpreal,0.0)
  168            elseif (arg.eq.'--magfield') then
  169               call getarg(k,arg)
  170               k=k+1
  171               j1=index(arg,':')
  172               read(arg(1:j1-1),*) magfieldstart
  173               arg(1:j1)=' '
  174               j2=index(arg,':')
  175               read(arg(j1+1:j2-1),*) magfieldstop
  176               read(arg(j2+1:),*) num_magfield_steps
  177            elseif (arg.eq.'--alpha') then
  178               call getarg(k,arg)
  179               k=k+1
  180               j1=index(arg,':')
  181               read(arg(1:j1-1),*) alphastart
  182               arg(1:j1)=' '
  183               j2=index(arg,':')
  184               read(arg(j1+1:j2-1),*) alphastop
  185               read(arg(j2+1:),*) num_alpha_steps
  186            elseif (arg.eq.'--beta') then
  187               call getarg(k,arg)
  188               k=k+1
  189               j1=index(arg,':')
  190               read(arg(1:j1-1),*) betastart
  191               arg(1:j1)=' '
  192               j2=index(arg,':')
  193               read(arg(j1+1:j2-1),*) betastop
  194               read(arg(j2+1:),*) num_beta_steps
  195            elseif (arg.eq.'--gamma') then
  196               call getarg(k,arg)
  197               k=k+1
  198               j1=index(arg,':')
  199               read(arg(1:j1-1),*) gammastart
  200               arg(1:j1)=' '
  201               j2=index(arg,':')
  202               read(arg(j1+1:j2-1),*) gammastop
  203               read(arg(j2+1:),*) num_gamma_steps
  204            elseif (arg.eq.'--delta') then
  205               call getarg(k,arg)
  206               k=k+1
  207               j1=index(arg,':')
  208               read(arg(1:j1-1),*) deltastart
  209               arg(1:j1)=' '
  210               j2=index(arg,':')
  211               read(arg(j1+1:j2-1),*) deltastop
  212               read(arg(j2+1:),*) num_delta_steps
  213            elseif (arg.eq.'--eta') then
  214               call getarg(k,arg)
  215               k=k+1
  216               j1=index(arg,':')
  217               read(arg(1:j1-1),*) etastart
  218               arg(1:j1)=' '
  219               j2=index(arg,':')
  220               read(arg(j1+1:j2-1),*) etastop
  221               read(arg(j2+1:),*) num_eta_steps
  222            elseif (arg.eq.'--kappa') then
  223               call getarg(k,arg)
  224               k=k+1
  225               j1=index(arg,':')
  226               read(arg(1:j1-1),*) kappastart
  227               arg(1:j1)=' '
  228               j2=index(arg,':')
  229               read(arg(j1+1:j2-1),*) kappastop
  230               read(arg(j2+1:),*) num_kappa_steps
  231            elseif (arg.eq.'--length') then
  232               call getarg(k,arg)
  233               k=k+1
  234               j1=index(arg,':')
  235               read(arg(1:j1-1),*) lengthstart
  236               arg(1:j1)=' '
  237               j2=index(arg,':')
  238               read(arg(j1+1:j2-1),*) lengthstop
  239               read(arg(j2+1:),*) num_length_steps
  240            elseif (arg.eq.'--zeta') then
  241               call getarg(k,arg)
  242               k=k+1
  243               j1=index(arg,':')
  244               read(arg(1:j1-1),*) zetastart
  245               arg(1:j1)=' '
  246               j2=index(arg,':')
  247               read(arg(j1+1:j2-1),*) zetastop
  248               read(arg(j2+1:),*) num_zeta_steps
  249            else
  250               write(6,20) 'Error: Unrecognized command line option: ',arg
  251    20         format('mopagras: ',a41,a26)
  252               stop
  253            endif
  254         enddo
  255   
  256   c     *********************************************************************
  257   c     Initialize the step lengths to be taken for each scan parameter.
  258   c     *********************************************************************
  259         if (num_magfield_steps.gt.1) then
  260            dmagfield=(magfieldstop-magfieldstart)/(num_magfield_steps-1)
  261            num_scan_params=num_scan_params+1
  262         else
  263            dzeta=(zetastop-zetastart)
  264         endif
  265         if (num_alpha_steps.gt.1) then
  266            dalpha=(alphastop-alphastart)/(num_alpha_steps-1)
  267            num_scan_params=num_scan_params+1
  268         else
  269            dalpha=(alphastop-alphastart)
  270         endif
  271         if (num_beta_steps.gt.1) then
  272            dbeta=(betastop-betastart)/(num_beta_steps-1)
  273            num_scan_params=num_scan_params+1
  274         else
  275            dbeta=(betastop-betastart)
  276         endif
  277         if (num_gamma_steps.gt.1) then
  278            dgamma=(gammastop-gammastart)/(num_gamma_steps-1)
  279            num_scan_params=num_scan_params+1
  280         else
  281            dgamma=(gammastop-gammastart)
  282         endif
  283         if (num_delta_steps.gt.1) then
  284            ddelta=(deltastop-deltastart)/(num_delta_steps-1)
  285            num_scan_params=num_scan_params+1
  286         else
  287            ddelta=(deltastop-deltastart)
  288         endif
  289         if (num_eta_steps.gt.1) then
  290            deta=(etastop-etastart)/(num_eta_steps-1)
  291            num_scan_params=num_scan_params+1
  292         else
  293            deta=(etastop-etastart)
  294         endif
  295         if (num_kappa_steps.gt.1) then
  296            dkappa=(kappastop-kappastart)/(num_kappa_steps-1)
  297            num_scan_params=num_scan_params+1
  298         else
  299            dkappa=(kappastop-kappastart)
  300         endif
  301         if (num_length_steps.gt.1) then
  302            dlength=(lengthstop-lengthstart)/(num_length_steps-1)
  303            num_scan_params=num_scan_params+1
  304         else
  305            dlength=(lengthstop-lengthstart)
  306         endif
  307         if (num_zeta_steps.gt.1) then
  308            dzeta=(zetastop-zetastart)/(num_zeta_steps-1)
  309            num_scan_params=num_scan_params+1
  310         else
  311            dzeta=(zetastop-zetastart)
  312         endif
  313         if (num_scan_params.gt.1) then
  314            write(6,20) 'error: more than one scan parameter specified!'
  315            write(6,20) '(please check parameter intervals for errors)'
  316            stop
  317         endif
  318   
  319   c     *********************************************************************
  320   c     If in verbose mode, show the parsed scan parameter interals at
  321   c     standard terminal output.
  322   c     *********************************************************************
  323         if (verbose.eqv.true) then
  324            write(6,30) '------------------------------------------------'
  325    30      format('------------',a,'------------')
  326            write(6,40) 'Summary of parameters parsed from command line'
  327    40      format(a)
  328            write(6,50) numargs
  329    50      format('number of arguments parsed: ',i2)
  330            write(6,60) 'Input signal field amplitude (lcp): ',a1iplus
  331            write(6,60) 'Input signal field amplitude (rcp): ',a1iminus
  332            write(6,60) 'Input idler field amplitude (lcp):  ',a2iplus
  333            write(6,60) 'Input idler field amplitude (rcp):  ',a2iminus
  334    60      format(3x,a36,e15.7,e15.7)
  335            write(6,70) 'Start value','Stop value','Increment','Scanning'
  336    70      format(10x,a16,a16,a16,a13)
  337            arg='(no scan)'
  338            if (num_magfield_steps.gt.1) arg='(scan)'
  339            write(6,80) 'magfield: ',magfieldstart,magfieldstop,
  340        +        dmagfield,arg
  341            arg='(no scan)'
  342            if (num_alpha_steps.gt.1) arg='(scan)'
  343            write(6,80) 'alpha:    ',alphastart,alphastop,dalpha,arg
  344            arg='(no scan)'
  345            if (num_beta_steps.gt.1) arg='(scan)'
  346            write(6,80) 'beta:     ',betastart,betastop,dbeta,arg
  347            arg='(no scan)'
  348            if (num_gamma_steps.gt.1) arg='(scan)'
  349            write(6,80) 'gamma:    ',gammastart,gammastop,dgamma,arg
  350            arg='(no scan)'
  351            if (num_delta_steps.gt.1) arg='(scan)'
  352            write(6,80) 'delta:    ',deltastart,deltastop,ddelta,arg
  353            arg='(no scan)'
  354            if (num_eta_steps.gt.1) arg='(scan)'
  355            write(6,80) 'eta:      ',etastart,etastop,deta,arg
  356            arg='(no scan)'
  357            if (num_kappa_steps.gt.1) arg='(scan)'
  358            write(6,80) 'kappa:    ',kappastart,kappastop,dkappa,arg
  359            arg='(no scan)'
  360            if (num_length_steps.gt.1) arg='(scan)'
  361            write(6,80) 'length:   ',lengthstart,lengthstop,dlength,arg
  362            arg='(no scan)'
  363            if (num_zeta_steps.gt.1) arg='(scan)'
  364            write(6,80) 'zeta:     ',zetastart,zetastop,dzeta,arg
  365    80      format(3x,a10,e15.7,'  ',e15.7,'  ',e15.7,' ',a9)
  366            write(6,30) '------------------------------------------------'
  367         endif
  368   
  369   c     *********************************************************************
  370   c     Open the file specified by the string ouputfilename for output.
  371   c     *********************************************************************
  372         open(oid_lcp,file=outfilename_lcp)
  373         open(oid_rcp,file=outfilename_rcp)
  374   
  375   c     *********************************************************************
  376   c     Enter the main loop of parameter scanning in the calculation of
  377   c     grating field envelopes.
  378   c     *********************************************************************
  379         do 100,i_alpha=1,num_alpha_steps
  380            alpha=alphastart+(i_alpha-1)*dalpha
  381            do 110,i_beta=1,num_beta_steps
  382               beta=betastart+(i_beta-1)*dbeta
  383               do 120,i_gamma=1,num_gamma_steps
  384                  gamma=gammastart+(i_gamma-1)*dgamma
  385                  do 130,i_delta=1,num_delta_steps
  386                     delta=deltastart+(i_delta-1)*ddelta
  387                     do 140,i_eta=1,num_eta_steps
  388                        eta=etastart+(i_eta-1)*deta
  389                        do 150,i_kappa=1,num_kappa_steps
  390                           kappa=kappastart+(i_kappa-1)*dkappa
  391                           do 160,i_length=1,num_length_steps
  392                              length=lengthstart+(i_length-1)*dlength
  393                              do 170,i_zeta=1,num_zeta_steps
  394                                 zeta=zetastart+(i_zeta-1)*dzeta
  395                                 do 180,i_magfield=1,num_magfield_steps
  396                                    magfield=magfieldstart
  397        +                                +(i_magfield-1)*dmagfield
  398                              if (num_magfield_steps.gt.1) then
  399                                 tmpgamma=magfield*gamma
  400                                 tmpdelta=magfield*delta
  401                              call mopagrat(alpha,beta,tmpgamma,
  402        +                       tmpdelta,eta,kappa,length,zeta,
  403        +                       a1iplus,a2iplus,a1iminus,a2iminus,
  404        +                       a1famplplus,a1bamplplus,a2famplplus,
  405        +                       a1famplminus,a1bamplminus,a2famplminus,
  406        +                       a1fplus,a1bplus,a2fplus,
  407        +                       a1fminus,a1bminus,a2fminus)
  408                              else
  409                              call mopagrat(alpha,beta,gamma,
  410        +                       delta,eta,kappa,length,zeta,
  411        +                       a1iplus,a2iplus,a1iminus,a2iminus,
  412        +                       a1famplplus,a1bamplplus,a2famplplus,
  413        +                       a1famplminus,a1bamplminus,a2famplminus,
  414        +                       a1fplus,a1bplus,a2fplus,
  415        +                       a1fminus,a1bminus,a2fminus)
  416                              endif
  417   
  418                              if (num_alpha_steps.gt.1) then
  419                                 write(oid_lcp,90) alpha,
  420        +                        a1famplplus,a1bamplplus,a2famplplus
  421                                 write(oid_rcp,90) alpha,
  422        +                        a1famplminus,a1bamplminus,a2famplminus
  423                              endif
  424                              if (num_beta_steps.gt.1) then
  425                                 write(oid_lcp,90) beta,
  426        +                        a1famplplus,a1bamplplus,a2famplplus
  427                                 write(oid_rcp,90) beta,
  428        +                        a1famplminus,a1bamplminus,a2famplminus
  429                              endif
  430                              if (num_magfield_steps.gt.1) then
  431                                 write(oid_lcp,90) magfield,
  432        +                        a1famplplus,a1bamplplus,a2famplplus
  433                                 write(oid_rcp,90) magfield,
  434        +                        a1famplminus,a1bamplminus,a2famplminus
  435                              else
  436                                 if (num_gamma_steps.gt.1) then
  437                                    write(oid_lcp,90) gamma,
  438        +                        a1famplplus,a1bamplplus,a2famplplus
  439                                    write(oid_rcp,90) gamma,
  440        +                        a1famplminus,a1bamplminus,a2famplminus
  441                                 endif
  442                                 if (num_delta_steps.gt.1) then
  443                                    write(oid_lcp,90) delta,
  444        +                        a1famplplus,a1bamplplus,a2famplplus
  445                                    write(oid_rcp,90) delta,
  446        +                        a1famplminus,a1bamplminus,a2famplminus
  447                                 endif
  448                              endif
  449                              if (num_eta_steps.gt.1) then
  450                                 write(oid_lcp,90) eta,
  451        +                        a1famplplus,a1bamplplus,a2famplplus
  452                                 write(oid_rcp,90) eta,
  453        +                        a1famplminus,a1bamplminus,a2famplminus
  454                              endif
  455                              if (num_kappa_steps.gt.1) then
  456                                 write(oid_lcp,90) kappa,
  457        +                        a1famplplus,a1bamplplus,a2famplplus
  458                                 write(oid_rcp,90) kappa,
  459        +                        a1famplminus,a1bamplminus,a2famplminus
  460                              endif
  461                              if (num_length_steps.gt.1) then
  462                                 write(oid_lcp,90) length,
  463        +                        a1famplplus,a1bamplplus,a2famplplus
  464                                 write(oid_rcp,90) length,
  465        +                        a1famplminus,a1bamplminus,a2famplminus
  466                              endif
  467                              if (num_zeta_steps.gt.1) then
  468                                 write(oid_lcp,90) zeta,
  469        +                        a1famplplus,a1bamplplus,a2famplplus
  470                                 write(oid_rcp,90) zeta,
  471        +                        a1famplminus,a1bamplminus,a2famplminus
  472                              endif
  473    180                          continue
  474    170                       continue
  475    160                    continue
  476    150                 continue
  477    140              continue
  478    130           continue
  479    120        continue
  480    110     continue
  481    100  continue
  482    90   format(1p,4(1x,e16.8))
  483   
  484   c     *********************************************************************
  485   c     Close the output files.
  486   c     *********************************************************************
  487         close(oid_lcp)
  488         close(oid_rcp)
  489   
  490         stop
  491         end
  492   

Return to previous page

Generated by ::viewsrc::

Last modified Wednesday 15 Feb 2023