      SUBROUTINE AHWFWC                                                 WFWC001

C          CALCULATES THE REFRACTIVE INDEX AND ITS GRADIENT USING THE   WFWC002

C          APPLETON-HARTREE FORMULA WITH FIELD, WITH COLLISIONS         WFWC003

      COMMON /CONST/ PI,PIT2,PID2,DEGS,RADIAN,K,C,LOGTEN                WFWC004

      COMMON /RIN/ MODRIN(3),COLL,FIELD,SPACE,KAY2,H,PHPT,PHPR,PHPTH,   WFWC005

     1             PHPPH,PHPOM,PHPKR,PHPKTH,PHPKPH,KPHPK,POLAR,LPOLAR   WFWC006

     2      ,SGN                                                        WFWC0065

      COMMON /XX/ MODX(2),X,PXPR,PXPTH,PXPPH,PXPT,HMAX                  WFWC007

      COMMON /YY/ MODY,Y,PYPR,PYPTH,PYPPH,YR,PYRPR,PYRPT,PYRPP,YTH,PYTPRWFWC008

     1            ,PYTPT,PYTPP,YPH,PYPPR,PYPPT,PYPPP                    WFWC009

      COMMON /ZZ/ MODZ,Z,PZPR,PZPTH,PZPPH                               WFWC010

      COMMON R,TH,PH,KR,KTH,KPH   /WW/ ID(10),W0,W(400)                 WFWC011

      COMMON /RK/ N,STEP,MODE,E1MAX,E1MIN,E2MAX,E2MIN,FACT,RSTART       WFWC012

      EQUIVALENCE (RAY,W(1)),(F,W(6))                                   WFWC013

      LOGICAL SPACE                                                     WFWC014

      REAL KR,KTH,KPH,K2                                                WFWC015

      COMPLEX N2,PNPR,PNPTH,PNPPH,PNPVR,PNPVTH,PNPVPH,NNP,PNPT,         WFWC016

     1        POLAR,LPOLAR,I,U,RAD,D,PNPPS,PNPX,PNPY,PNPZ,UX,UX2,D2,    WFWC017

     2        KAY2,H,PHPT,PHPR,PHPTH,PHPPH,PHPOM,PHPKR,PHPKTH,PHPKPH,   WFWC018

     3        KPHPK                                                     WFWC019

      DATA (MODRIN=8HAPPLETON,8H-HARTREE,8H FORMULA),(COLL=1.),         WFWC020

     1     (FIELD=1.),                                                  WFWC021

     2     (X=0.),(PXPR=0.),(PXPTH=0.),(PXPPH=0.),(PXPT=0.),            WFWC022

     3     (Y=0.),(PYPR=0.),(PYPTH=0.),(PYPPH=0.),(YR=0.),(PYRPR=0.),   WFWC023

     4     (PYRPT=0.),(PYRPP=0.),(YTH=0.),(PYTPR=0.),(PYTPT=0.),        WFWC024

     5     (PYTPP=0.),(YPH=0.),(PYPPR=0.),(PYPPT=0.),(PYPPP=0.)         WFWC025

     6    ,(Z=0.),(PZPR=0.),(PZPTH=0.),(PZPPH=0.),                      WFWC026

     7     (I=(0.,1.)),(ABSLIM=1.E-5)                                   WFWC027

      ENTRY RINDEX                                                      WFWC028

      OM=PIT2*1.E6*F                                                    WFWC029

      C2=C*C                                                            WFWC030

      K2=KR*KR+KTH*KTH+KPH*KPH                                          WFWC031

      OM2=OM*OM                                                         WFWC032

      VR =C/OM*KR                                                       WFWC033

      VTH=C/OM*KTH                                                      WFWC034

      VPH=C/OM*KPH                                                      WFWC035

      CALL ELECTX                                                       WFWC036

      CALL MAGY                                                         WFWC037

      V2=VR**2+VTH**2+VPH**2                                            WFWC038

      VDOTY=VR*YR+VTH*YTH+VPH*YPH                                       WFWC039

      YLV=VDOTY/V2                                                      WFWC040

      YL2=VDOTY**2/V2                                                   WFWC041

      YT2=Y**2-YL2                                                      WFWC042

      YT4=YT2*YT2                                                       WFWC043

      CALL COLFRZ                                                       WFWC044

      U=CMPLX(1.,-Z)                                                    WFWC045

      UX=U-X                                                            WFWC046

      UX2=UX*UX                                                         WFWC047

      RAD=RAY*CSQRT(YT4+4.*YL2*UX2)                                     WFWC048

      D=2.*U*UX-YT2+RAD                                                 WFWC049

      D2=D*D                                                            WFWC050

      N2=1.-2.*X*UX/D                                                   WFWC051

      PNPPS=2.*X*UX*(-1.+(YT2-2.*UX2)/RAD)/D2                           WFWC052

      PPSPR =YL2/Y*PYPR -(VR*PYRPR+VTH*PYTPR+VPH*PYPPR)*YLV             WFWC053

      PPSPTH=YL2/Y*PYPTH-(VR*PYRPT+VTH*PYTPT+VPH*PYPPT)*YLV             WFWC054

      PPSPPH=YL2/Y*PYPPH-(VR*PYRPP+VTH*PYTPP+VPH*PYPPP)*YLV             WFWC055

      PNPX=-(2.*U*UX2-YT2*(U-2.*X)+(YT4*(U-2.*X)+4.*YL2*UX*UX2)/RAD)/D2 WFWC056

      PNPY=2.*X*UX*(-YT2+(YT4+2.*YL2*UX2)/RAD)/(D2*Y)                   WFWC057

      PNPZ=I*X*(-2.*UX2-YT2+YT4/RAD)/D2                                 WFWC058

      PNPR =PNPX*PXPR +PNPY*PYPR +PNPZ*PZPR +PNPPS*PPSPR                WFWC059

      PNPTH=PNPX*PXPTH+PNPY*PYPTH+PNPZ*PZPTH+PNPPS*PPSPTH               WFWC060

      PNPPH=PNPX*PXPPH+PNPY*PYPPH+PNPZ*PZPPH+PNPPS*PPSPPH               WFWC061

      PNPVR =PNPPS*(VR *YL2/V2-YLV*YR )                                 WFWC062

      PNPVTH=PNPPS*(VTH*YL2/V2-YLV*YTH)                                 WFWC063

      PNPVPH=PNPPS*(VPH*YL2/V2-YLV*YPH)                                 WFWC064

      NNP=N2-(2.*X*PNPX+Y*PNPY+Z*PNPZ)                                  WFWC065

      PNPT=PNPX*PXPT                                                    WFWC066

      SPACE=REAL(N2).EQ.1..AND.ABS(AIMAG(N2)).LT.ABSLIM                 WFWC067

      POLAR=-I*SQRT(V2)*(-YT2+RAD)/(2.*VDOTY*UX)                        WFWC068

      GAM=(-YT2+RAD)/(2.*UX)                                            WFWC069

      LPOLAR=I*X*SQRT(YT2)/(UX*(U+GAM))                                 WFWC070

      KAY2=OM2/C2*N2                                                    WFWC071

      IF(RSTART.EQ.0.) GO TO 1                                          WFWC072

      SCALE=SQRT(REAL(KAY2)/K2)                                         WFWC073

      KR =SCALE*KR                                                      WFWC074

      KTH=SCALE*KTH                                                     WFWC075

      KPH=SCALE*KPH                                                     WFWC076

  1   CONTINUE                                                          WFWC077

C********* CALCULATES A HAMILTONIAN H                                   WFWC078

      H=.5*(C2*K2/OM2-N2)                                               WFWC079

C********* AND ITS PARTIAL DERIVATIVES WITH RESPECT TO                  WFWC080

C********* TIME, R, THETA, PHI, OMEGA, KR, KTHETA, AND KPHI.            WFWC081

      PHPT =-PNPT                                                       WFWC082

      PHPR =-PNPR                                                       WFWC083

      PHPTH=-PNPTH                                                      WFWC084

      PHPPH=-PNPPH                                                      WFWC085

      PHPOM=-NNP/OM                                                     WFWC086

      PHPKR =C2/OM2*KR -C/OM*PNPVR                                      WFWC087

      PHPKTH=C2/OM2*KTH-C/OM*PNPVTH                                     WFWC088

      PHPKPH=C2/OM2*KPH-C/OM*PNPVPH                                     WFWC089

      KPHPK=N2                                                          WFWC090

      RETURN                                                            WFWC091

         END                                                            WFWC092

      SUBROUTINE AHWFNC                                                 WFNC001

C          CALCULATES THE REFRACTIVE INDEX AND ITS GRADIENT USING THE   WFNC002

C          APPLETON-HARTREE FORMULA WITH FIELD, NO COLLISIONS           WFNC003

      COMMON /CONST/ PI,PIT2,PID2,DEGS,RADIAN,K,C,LOGTEN                WFNC004

      COMMON /RIN/ MODRIN(3),COLL,FIELD,SPACE,KAY2,KAY2I,               WFNC005

     1              H,HI,PHPT,PHPTI,PHPR,PHPRI,PHPTH,PHPTHI,PHPPH,PHPPHIWFNC006

     2            ,PHPOM,PHPOMI,PHPKR,PHPKRI,PHPKTH,PHPKTI,PHPKPH,PHPKPIWFNC007

     3            ,KPHPK,KPHPKI,POLAR,POLARI,LPOLAR,LPOLRI,SGN          WFNC008

      COMMON /XX/ MODX(2),X,PXPR,PXPTH,PXPPH,PXPT,HMAX                  WFNC009

      COMMON /YY/ MODY,Y,PYPR,PYPTH,PYPPH,YR,PYRPR,PYRPT,PYRPP,YTH,PYTPRWFNC010

     1           ,PYTPT,PYTPP,YPH,PYPPR,PYPPT,PYPPP                     WFNC011

      COMMON /ZZ/ MODZ,Z(4)                                             WFNC012

      COMMON /RK/ N,STEP,MODE,E1MAX,E1MIN,E2MAX,E2MIN,FACT,RSTART       WFNC013

      COMMON R,TH,PH,KR,KTH,KPH   /WW/ ID(10),W0,W(400)                 WFNC014

      EQUIVALENCE (RAY,W(1)),(F,W(6))                                   WFNC015

      LOGICAL SPACE                                                     WFNC016

      REAL KR,KTH,KPH,K2,KPHPK,KPHPKI,KAY2,KAY2I,N2,NNP,LPOLAR,LPOLRI   WFNC017

      DATA (MODRIN=8HAPPLETON,8H-HARTREE,8H FORMULA),(COLL=0.),         WFNC018

     1     (FIELD=1.),(KAY2I=0.),(HI=0.),(PHPTI=0.),(PHPRI=0.),         WFNC019

     2     (PHPTHI=0.),(PHPPHI=0.),(PHPOMI=0.),(PHPKRI=0.),(PHPKTI=0.), WFNC020

     3     (PHPKPI=0.),(KPHPKI=0.),(POLAR=0.),(LPOLAR=0.),              WFNC021

     4     (X=0.),(PXPR=0.),(PXPTH=0.),(PXPPH=0.),(PXPT=0.),            WFNC022

     5     (Y=0.),(PYPR=0.),(PYPTH=0.),(PYPPH=0.),(YR=0.),(PYRPR=0.),   WFNC023

     6     (PYRPT=0.),(PYRPP=0.),(YTH=0.),(PYTPR=0.),(PYTPT=0.),        WFNC024

     7     (PYTPP=0.),(YPH=0.),(PYPPR=0.),(PYPPT=0.),(PYPPP=0.),        WFNC025

     8     (MODZ=1H ),(U=1.)                                            WFNC026

      ENTRY RINDEX                                                      WFNC027

      OM=PIT2*1.E6*F                                                    WFNC028

      C2=C*C                                                            WFNC029

      K2=KR*KR+KTH*KTH+KPH*KPH                                          WFNC030

      OM2=OM*OM                                                         WFNC031

      VR =C/OM*KR                                                       WFNC032

      VTH=C/OM*KTH                                                      WFNC033

      VPH=C/OM*KPH                                                      WFNC034

      CALL ELECTX                                                       WFNC035

      CALL MAGY                                                         WFNC036

      V2=VR**2+VTH**2+VPH**2                                            WFNC037

      VDOTY=VR*YR+VTH*YTH+VPH*YPH                                       WFNC038

      YLV=VDOTY/V2                                                      WFNC039

      YL2=VDOTY**2/V2                                                   WFNC040

      YT2=Y**2-YL2                                                      WFNC041

      YT4=YT2*YT2                                                       WFNC042

      UX=U-X                                                            WFNC043

      UX2=UX*UX                                                         WFNC044

      RAD=RAY*SQRT(YT4+4.*YL2*UX2)                                      WFNC045

      D=2.*UX-YT2+RAD                                                   WFNC046

      D2=D*D                                                            WFNC047

      N2=1.-2.*X*UX/D                                                   WFNC048

      PNPPS=2.*X*UX*(-1.+(YT2-2.*UX2)/RAD)/D2                           WFNC049

      PPSPR= YL2/Y*PYPR -(VR*PYRPR+VTH*PYTPR+VPH*PYPPR)*YLV             WFNC050

      PPSPTH=YL2/Y*PYPTH-(VR*PYRPT+VTH*PYTPT+VPH*PYPPT)*YLV             WFNC051

      PPSPPH=YL2/Y*PYPPH-(VR*PYRPP+VTH*PYTPP+VPH*PYPPP)*YLV             WFNC052

      PNPX=-(2.*UX2-YT2*(U-2.*X)+(YT4*(U-2.*X)+4.*YL2*UX*UX2)/RAD)/D2   WFNC053

      PNPY=2.*X*UX*(-YT2+(YT4+2.*YL2*UX2)/RAD)/(D2*Y)                   WFNC054

      NNP=N2-(2.*X*PNPX+Y*PNPY)                                         WFNC055

      PNPR =PNPX*PXPR +PNPY*PYPR +PNPPS*PPSPR                           WFNC056

      PNPTH=PNPX*PXPTH+PNPY*PYPTH+PNPPS*PPSPTH                          WFNC057

      PNPPH=PNPX*PXPPH+PNPY*PYPPH+PNPPS*PPSPPH                          WFNC058

      PNPVR =PNPPS*(VR *YL2-VDOTY*YR )/V2                               WFNC059

      PNPVTH=PNPPS*(VTH*YL2-VDOTY*YTH)/V2                               WFNC060

      PNPVPH=PNPPS*(VPH*YL2-VDOTY*YPH)/V2                               WFNC061

      PNPT=PNPX*PXPT                                                    WFNC062

      SPACE=N2.EQ.1.                                                    WFNC063

      POLARI=SQRT(V2)*(YT2-RAD)/(2.*VDOTY*UX)                           WFNC064

      GAM=(-YT2+RAD)/(2.*UX)                                            WFNC065

      LPOLRI=X*SQRT(YT2)/(UX*(U+GAM))                                   WFNC066

      KAY2=OM2/C2*N2                                                    WFNC067

      IF(RSTART.EQ.0.) GO TO 1                                          WFNC068

      SCALE=SQRT(KAY2/K2)                                               WFNC069

      KR =SCALE*KR                                                      WFNC070

      KTH=SCALE*KTH                                                     WFNC071

      KPH=SCALE*KPH                                                     WFNC072

  1   CONTINUE                                                          WFNC073

C********* CALCULATES A HAMILTONIAN H                                   WFNC074

      H=.5*(C2*K2/OM2-N2)                                               WFNC075

C********* AND ITS PARTIAL DERIVATIVES WITH RESPECT TO                  WFNC076

C********* TIME, R, THETA, PHI, OMEGA, KR, KTHETA, AND KPHI.            WFNC077

      PHPT =-PNPT                                                       WFNC078

      PHPR =-PNPR                                                       WFNC079

      PHPTH=-PNPTH                                                      WFNC080

      PHPPH=-PNPPH                                                      WFNC081

      PHPOM=-NNP/OM                                                     WFNC082

      PHPKR =C2/OM2*KR -C/OM*PNPVR                                      WFNC083

      PHPKTH=C2/OM2*KTH-C/OM*PNPVTH                                     WFNC084

      PHPKPH=C2/OM2*KPH-C/OM*PNPVPH                                     WFNC085

      KPHPK=N2                                                          WFNC086

      RETURN                                                            WFNC087

         END                                                            WFNC088-

      SUBROUTINE AHNFWC                                                 NFWC010

C          CALCULATES THE REFRACTIVE INDEX AND ITS GRADIENT USING THE   NFWC011

C          APPLETON-HARTREE FORMULA -- NO FIELD, WITH COLLISIONS        NFWC012

      COMMON /CONST/ PI,PIT2,PID2,DEGS,RADIAN,K,C,LOGTEN                NFWC013

      COMMON /RIN/ MODRIN(3),COLL,FIELD,SPACE,KAY2,H,PHPT,PHPR,PHPTH,   NFWC014

     1             PHPPH,PHPOM,PHPKR,PHPKTH,PHPKPH,KPHPK,POLAR,LPOLAR,  NFWC015

     2             SGN                                                  NFWC016

      COMMON /XX/ MODX(2),X,PXPR,PXPTH,PXPPH,PXPT,HMAX                  NFWC017

      COMMON /YY/ MODY,Y(16)                                            NFWC018

      COMMON /ZZ/ MODZ,Z,PZPR,PZPTH,PZPPH                               NFWC019

      COMMON /RK/ N,STEP,MODE,E1MAX,E1MIN,E2MAX,E2MIN,FACT,RSTART       NFWC020

      COMMON R,TH,PH,KR,KTH,KPH   /WW/ ID(10),W0,W(400)                 NFWC021

      EQUIVALENCE (RAY,W(1)),(F,W(6))                                   NFWC022

      LOGICAL SPACE                                                     NFWC023

      REAL KR,KTH,KPH,K2                                                NFWC024

      COMPLEX KAY2,H,PHPT,PHPR,PHPTH,PHPPH,PHPOM,PHPKR,PHPKTH,PHPKPH,   NFWC025

     1        KPHPK,POLAR,LPOLAR,U,I,PNPX,PNPZ,                         NFWC026

     2        N2,PNPR,PNPTH,PNPPH,PNPVR,PNPVTH,PNPVPH,NNP,PNPT          NFWC027

      DATA (MODRIN=8HAPPLETON,8H-HARTREE,8H FORMULA),(COLL=1.),         NFWC028

     1     (FIELD=0.),(POLAR=(0.,1.)),(LPOLAR=(0.,0.)),                 NFWC029

     2     (X=0.),(PXPR=0.),(PXPTH=0.),(PXPPH=0.),(PXPT=0.),            NFWC030

     3     (MODY=1H ),                                                  NFWC031

     4     (Z=0.),(PZPR=0.),(PZPTH=0.),(PZPPH=0.),                      NFWC032

     5     (I=(0.,1.)),(ABSLIM=1.E-5),(PNPVR=0.),(PNPVTH=0.),(PNPVPH=0.)NFWC033

      ENTRY RINDEX                                                      NFWC034

      OM=PIT2*1.E6*F                                                    NFWC035

      C2=C*C                                                            NFWC036

      K2=KR*KR+KTH*KTH+KPH*KPH                                          NFWC037

      OM2=OM*OM                                                         NFWC038

      VR =C/OM*KR                                                       NFWC039

      VTH=C/OM*KTH                                                      NFWC040

      VPH=C/OM*KPH                                                      NFWC041

      CALL ELECTX                                                       NFWC042

      CALL COLFRZ                                                       NFWC043

      U=1.-I*Z                                                          NFWC044

      N2=1.-X/U                                                         NFWC045

      PNPX=-1./(2.*U)                                                   NFWC046

      PNPZ=-I*X/(2.*U**2)                                               NFWC047

      NNP=N2-(2.*X*PNPX+Z*PNPZ)                                         NFWC048

      PNPR =PNPX*PXPR +PNPZ*PZPR                                        NFWC049

      PNPTH=PNPX*PXPTH+PNPZ*PZPTH                                       NFWC050

      PNPPH=PNPX*PXPPH+PNPZ*PZPPH                                       NFWC051

      PNPT=PNPX*PXPT                                                    NFWC052

      SPACE=REAL(N2).EQ.1..AND.ABS(AIMAG(N2)).LT.ABSLIM                 NFWC053

      KAY2=OM2/C2*N2                                                    NFWC054

      IF(RSTART.EQ.0.) GO TO 1                                          NFWC055

      SCALE=SQRT(REAL(KAY2)/K2)                                         NFWC056

      KR =SCALE*KR                                                      NFWC057

      KTH=SCALE*KTH                                                     NFWC058

      KPH=SCALE*KPH                                                     NFWC059

  1   CONTINUE                                                          NFWC060

C********* CALCULATES A HAMILTONIAN H                                   NFWC061

      H=.5*(C2*K2/OM2-N2)                                               NFWC062

C********* AND ITS PARTIAL DERIVATIVES WITH RESPECT TO                  NFWC063

C********* TIME, R, THETA, PHI, OMEGA, KR, KTHETA, AND KPHI.            NFWC064

      PHPT =-PNPT                                                       NFWC065

      PHPR =-PNPR                                                       NFWC066

      PHPTH=-PNPTH                                                      NFWC067

      PHPPH=-PNPPH                                                      NFWC068

      PHPOM=-NNP/OM                                                     NFWC069

      PHPKR =C2/OM2*KR                                                  NFWC070

      PHPKTH=C2/OM2*KTH                                                 NFWC071

      PHPKPH=C2/OM2*KPH                                                 NFWC072

      KPHPK=N2                                                          NFWC073

      RETURN                                                            NFWC074

         END                                                            NFWC075-

      SUBROUTINE AHNFNC                                                 NFNC001

C          CALCULATES THE REFRACTIVE INDEX AND ITS GRADIENT USING THE   NFNC002

C          APPLETON-HARTREE FORMULA -- NO FIELD, NO COLLISIONS          NFNC003

      COMMON /CONST/ PI,PIT2,PID2,DEGS,RADIAN,K,C,LOGTEN                NFNC004

      COMMON /RIN/ MODRIN(3),COLL,FIELD,SPACE,KAY2,KAY2I,               NFNC005

     1             H,HI,PHPT,PHPTI,PHPR,PHPRI,PHPTH,PHPTHI,PHPPH,PHPPHI,NFNC006

     2             PHPOM,PHPOMI,PHPKR,PHPKRI,PHPKTH,PHPKTI,PHPKPH,PHPKPINFNC007

     3            ,KPHPK,KPHPKI,POLAR,POLARI,LPOLAR,LPOLRI,SGN          NFNC008

      COMMON /XX/ MODX(2),X,PXPR,PXPTH,PXPPH,PXPT,HMAX                  NFNC009

      COMMON /YY/ MODY,Y(16)  /ZZ/ MODZ,Z(4)                            NFNC010

      COMMON /RK/ N,STEP,MODE,E1MAX,E1MIN,E2MAX,E2MIN,FACT,RSTART       NFNC011

      COMMON R,TH,PH,KR,KTH,KPH   /WW/ ID(10),W0,W(400)                 NFNC012

      EQUIVALENCE (RAY,W(1)),(F,W(6))                                   NFNC013

      LOGICAL SPACE                                                     NFNC014

      REAL N2,NNP,KR,KTH,KPH,K2,KPHPK,KPHPKI,KAY2,KAY2I,LPOLAR,LPOLRI   NFNC015

      DATA (MODRIN=8HAPPLETON,8H-HARTREE,8H FORMULA),(COLL=0.),         NFNC016

     1     (FIELD=0.),(KAY2I=0.),(HI=0.),(PHPTI=0.),(PHPRI=0.),         NFNC017

     2     (PHPTHI=0.),(PHPPHI=0.),(PHPOMI=0.),(PHPKRI=0.),(PHPKTI=0.), NFNC018

     3     (PHPKPI=0.),(KPHPKI=0.),(POLAR=0.),(POLARI=1.),(LPOLAR=0.),  NFNC019

     4     (LPOLRI=1.),                                                 NFNC020

     4     (X=0.),(PXPR=0.),(PXPTH=0.),(PXPPH=0.),(PXPT=0.),            NFNC021

     5     (MODY=1H ),(MODZ=1H ),                                       NFNC022

     6     (NNP=1.),(PNPX=-0.5),(PNPVR=0.),(PNPVTH=0.),(PNPVPH=0.)      NFNC023

      ENTRY RINDEX                                                      NFNC024

      OM=PIT2*1.E6*F                                                    NFNC025

      C2=C*C                                                            NFNC026

      K2=KR*KR+KTH*KTH+KPH*KPH                                          NFNC027

      OM2=OM*OM                                                         NFNC028

      VR =C/OM*KR                                                       NFNC029

      VTH=C/OM*KTH                                                      NFNC030

      VPH=C/OM*KPH                                                      NFNC031

      CALL ELECTX                                                       NFNC032

      PNPR =PNPX*PXPR                                                   NFNC033

      PNPTH=PNPX*PXPTH                                                  NFNC034

      PNPPH=PNPX*PXPPH                                                  NFNC035

      PNPT=PNPX*PXPT                                                    NFNC036

      N2=1.-X                                                           NFNC037

      SPACE=N2.EQ.1.                                                    NFNC038

      KAY2=OM2/C2*N2                                                    NFNC039

      IF(RSTART.EQ.0.) GO TO 1                                          NFNC040

      SCALE=SQRT(KAY2/K2)                                               NFNC041

      KR =SCALE*KR                                                      NFNC042

      KTH=SCALE*KTH                                                     NFNC043

      KPH=SCALE*KPH                                                     NFNC044

  1   CONTINUE                                                          NFNC045

C********* CALCULATES A HAMILTONIAN H                                   NFNC046

      H=.5*(C2*K2/OM2-N2)                                               NFNC047

C********* AND ITS PARTIAL DERIVATIVES WITH RESPECT TO                  NFNC048

C********* TIME, R, THETA, PHI, OMEGA, KR, KTHETA, AND KPHI.            NFNC049

      PHPT =-PNPT                                                       NFNC050

      PHPR =-PNPR                                                       NFNC051

      PHPTH=-PNPTH                                                      NFNC052

      PHPPH=-PNPPH                                                      NFNC053

      PHPOM=-NNP/OM                                                     NFNC054

      PHPKR =C2/OM2*KR                                                  NFNC055

      PHPKTH=C2/OM2*KTH                                                 NFNC056

      PHPKPH=C2/OM2*KPH                                                 NFNC057

      KPHPK=N2                                                          NFNC058

      RETURN                                                            NFNC059

         END                                                            NFNC060

      SUBROUTINE BQWFWC                                                 BQWC001

C********* CALCULATES A HAMILTONIAN H                                   BQWC002

C********* (= BOOKER QUARTIC FOR VERTICAL INCIDENCE, S=0, C=1)          BQWC003

C********* AND ITS PARTIAL DERIVATIVES WITH RESPECT TO                  BQWC004

C********* TIME, R, THETA, PHI, OMEGA, KR, KTHETA, AND KPHI.            BQWC005

C********* WITH FIELD, WITH COLLISIONS                                  BQWC006

      COMMON /CONST/ PI,PIT2,PID2,DEGS,RADIAN,K,C,LOGTEN                BQWC007

      COMMON /RIN/ MODRIN(3),COLL,FIELD,SPACE,KAY2,H,PHPT,PHPR,PHPTH,   BQWC008

     1             PHPPH,PHPOM,PHPKR,PHPKTH,PHPKPH,KPHPK,POLAR,LPOLAR,  BQWC009

     2             SGN                                                  BQWC010

      COMMON /XX/ MODX(2),X,PXPR,PXPTH,PXPPH,PXPT,HMAX                  BQWC011

      COMMON /YY/ MODY,Y,PYPR,PYPTH,PYPPH,YR,PYRPR,PYRPT,PYRPP,YTH,PYTPRBQWC012

     1           ,PYTPT,PYTPP,YPH,PYPPR,PYPPT,PYPPP                     BQWC013

      COMMON /ZZ/ MODZ,Z,PZPR,PZPTH,PZPPH                               BQWC014

      COMMON R,TH,PH,KR,KTH,KPH   /WW/ ID(10),W0,W(400)                 BQWC015

      COMMON /RK/ N,STEP,MODE,E1MAX,E1MIN,E2MAX,E2MIN,FACT,RSTART       BQWC016

      COMMON /FLG/ NTYP,NEWWR,NEWWP,PENET,LINES,IHOP,HPUNCH             BQWC017

      EQUIVALENCE (RAY,W(1)),(F,W(6))                                   BQWC018

      LOGICAL SPACE                                                     BQWC019

      REAL KR,KTH,KPH,K2,KDOTY,K4,KDOTY2                                BQWC020

      COMPLEX KAY2,H,PHPT,PHPR,PHPTH,PHPPH,PHPOM,PHPKR,PHPKTH,PHPKPH,   BQWC021

     1        POLAR,LPOLAR,I,U,RAD,D,PNPPS,PNPX,PNPY,PNPZ,UX,UX2,D2,    BQWC022

     2        KPHPK,U2,A,B,ALPHA,BETA,GAMMA,PHPX,PHPY2,PHPK2,PHPU,PHPZ, BQWC023

     3        N2,PNPR,PNPTH,PNPPH,PNPVR,PNPVTH,PNPVPH,NNP,PNPT          BQWC024

      DATA (MODRIN=8HBOOKER Q,8HUARTIC, ,8HS=0, C=1),(COLL=1.),         BQWC025

     1     (FIELD=1.),                                                  BQWC026

     2     (X=0.),(PXPR=0.),(PXPTH=0.),(PXPPH=0.),(PXPT=0.),            BQWC027

     3     (Y=0.),(PYPR=0.),(PYPTH=0.),(PYPPH=0.),(YR=0.),(PYRPR=0.),   BQWC028

     4     (PYRPT=0.),(PYRPP=0.),(YTH=0.),(PYTPR=0.),(PYTPT=0.),        BQWC029

     5     (PYTPP=0.),(YPH=0.),(PYPPR=0.),(PYPPT=0.),(PYPPP=0.)         BQWC030

     6    ,(Z=0.),(PZPR=0.),(PZPTH=0.),(PZPPH=0.),                      BQWC031

     7     (I=(0.,1.)),(ABSLIM=1.E-5),(SGN=1.)                          BQWC032

      ENTRY RINDEX                                                      BQWC033

      OM=PIT2*1.E6*F                                                    BQWC034

      C2=C*C                                                            BQWC035

      K2=KR*KR+KTH*KTH+KPH*KPH                                          BQWC036

      OM2=OM*OM                                                         BQWC037

      CALL ELECTX                                                       BQWC038

      IF(X.LT..1) GO TO 2                                               BQWC039

      K4=K2*K2                                                          BQWC040

      OM4=OM2*OM2                                                       BQWC041

      C4=C2*C2                                                          BQWC042

      CALL MAGY                                                         BQWC043

      Y2=Y*Y                                                            BQWC044

      KDOTY=KR*YR+KTH*YTH+KPH*YPH                                       BQWC045

      KDOTY2=KDOTY*KDOTY                                                BQWC046

      CALL COLFRZ                                                       BQWC047

      U=CMPLX(1.,-Z)                                                    BQWC048

      U2=U*U                                                            BQWC049

      UX=U-X                                                            BQWC050

      UX2=UX*UX                                                         BQWC051

      A=UX*U2-U*Y2                                                      BQWC052

      B=-2.*U*UX2+Y2*(2.*U-X)                                           BQWC053

      ALPHA=A*C4*K4+X*KDOTY2*C4*K2                                      BQWC054

      BETA=B*C2*OM2*K2-X*KDOTY2*C2*OM2                                  BQWC055

      GAMMA=(UX2-Y2)*UX*OM4                                             BQWC056

      H=ALPHA+BETA+GAMMA                                                BQWC057

      PHPX=-U2*C4*K4+KDOTY2*C4*K2+(4.*U*UX-Y2)*C2*OM2*K2-KDOTY2*C2*OM2+ BQWC058

     1 (-3.*UX2+Y2)*OM4                                                 BQWC059

      PHPY2=-U*C4*K4+(2.*U-X)*C2*OM2*K2-UX*OM4                          BQWC060

      PHPKY2   =X*C2*(C2*K2-OM2)                                        BQWC061

      PHPU=(2.*U*UX+U2-Y2)*C4*K4+2.*(Y2-UX2-2.*U*UX)*C2*K2*OM2+(3.*UX2  BQWC062

     1 -Y2)*OM4                                                         BQWC063

      PHPZ=-I*PHPU                                                      BQWC064

      PHPK2=2.*A*C4*K2+X*KDOTY2*C4+B*C2*OM2                             BQWC065

      PHPT=PHPX*PXPT                                                    BQWC066

      PHPR =PHPX*PXPR +PHPY2*2.*Y*PYPR +PHPKY2  *2.*KDOTY*              BQWC067

     1 (KR*PYRPR+KTH*PYTPR+KPH*PYPPR) +PHPZ*PZPR                        BQWC068

      PHPTH=PHPX*PXPTH+PHPY2*2.*Y*PYPTH+PHPKY2  *2.*KDOTY*              BQWC069

     1 (KR*PYRPT+KTH*PYTPT+KPH*PYPPT) +PHPZ*PZPTH                       BQWC070

      PHPPH=PHPX*PXPPH+PHPY2*2.*Y*PYPPH+PHPKY2  *2.*KDOTY*              BQWC071

     1 (KR*PYRPP+KTH*PYTPP+KPH*PYPPP) +PHPZ*PZPPH                       BQWC072

      PHPOM=(2.*BETA+4.*GAMMA)/OM                                       BQWC073

     1 -2.*PHPX*X/OM-2.*PHPY2*Y2/OM-2.*PHPKY2  *KDOTY2/OM -PHPZ*Z/OM    BQWC074

      PHPKR= 2.*PHPK2*KR +2.*KDOTY*PHPKY2  *YR                          BQWC075

      PHPKTH=2.*PHPK2*KTH+2.*KDOTY*PHPKY2  *YTH                         BQWC076

      PHPKPH=2.*PHPK2*KPH+2.*KDOTY*PHPKY2  *YPH                         BQWC077

      KAY2=K2*(-BETA+SGN*RAY*CSQRT(BETA**2-4.*ALPHA*GAMMA))/(2.*ALPHA)  BQWC078

C                                                                       BQWC079

      IF(RSTART.EQ.0.) GO TO 1                                          BQWC080

      SCALE=SQRT((-REAL(BETA)+SGN*RAY*SQRT(REAL(BETA)**2                BQWC081

     1 -4.*REAL(ALPHA)*REAL(GAMMA)))/(2.*REAL(ALPHA)))                  BQWC082

      KR =SCALE*KR                                                      BQWC083

      KTH=SCALE*KTH                                                     BQWC084

      KPH=SCALE*KPH                                                     BQWC085

  1   CONTINUE                                                          BQWC086

C********* THE FOLLOWING 3 CARDS USED FOR RAY TRACING IN COMPLEX SPACE  BQWC087

C     IF(CABS((-BETA-SGN*RAY*CSQRT(BETA**2-4.*ALPHA*GAMMA))/ALPHA-2.).  BQWC088

C    1LT.CABS((-BETA+SGN*RAY*CSQRT(BETA**2-4.*ALPHA*GAMMA))/ALPHA-2.)   BQWC089

C    2 .AND.RSTART.EQ.0.) SGN=-SGN                                      BQWC090

      KPHPK=4.*ALPHA+2.*BETA                                            BQWC091

      SPACE=CABS(C2*KAY2/OM2-1.).LT.ABSLIM                              BQWC092

      POLAR =SQRT(K2)*(U+X*OM2/(C2*KAY2-OM2))/KDOTY*I                   BQWC093

      LPOLR = SQRT(Y2-KDOTY2/K2)/UX*(1.-C2*KAY2/OM2)*I                  BQWC094

      RETURN                                                            BQWC095

C          CALCULATES THE REFRACTIVE INDEX AND ITS GRADIENT USING THE   BQWC096

C          APPLETON-HARTREE FORMULA WITH FIELD, WITH COLLISIONS         BQWC097

  2   CONTINUE                                                          BQWC098

      VR =C/OM*KR                                                       BQWC099

      VTH=C/OM*KTH                                                      BQWC100

      VPH=C/OM*KPH                                                      BQWC101

      CALL MAGY                                                         BQWC102

      V2=VR**2+VTH**2+VPH**2                                            BQWC103

      VDOTY=VR*YR+VTH*YTH+VPH*YPH                                       BQWC104

      YLV=VDOTY/V2                                                      BQWC105

      YL2=VDOTY**2/V2                                                   BQWC106

      YT2=Y**2-YL2                                                      BQWC107

      YT4=YT2*YT2                                                       BQWC108

      CALL COLFRZ                                                       BQWC109

      U=CMPLX(1.,-Z)                                                    BQWC110

      UX=U-X                                                            BQWC111

      UX2=UX*UX                                                         BQWC112

      RAD=SGN*RAY*CSQRT(YT4+4.*YL2*UX2)                                 BQWC113

      D=2.*U*UX-YT2+RAD                                                 BQWC114

      D2=D*D                                                            BQWC115

      N2=1.-2.*X*UX/D                                                   BQWC116

      PNPPS=2.*X*UX*(-1.+(YT2-2.*UX2)/RAD)/D2                           BQWC117

      PPSPR= YL2/Y*PYPR -(VR*PYRPR+VTH*PYTPR+VPH*PYPPR)*YLV             BQWC118

      PPSPTH=YL2/Y*PYPTH-(VR*PYRPT+VTH*PYTPT+VPH*PYPPT)*YLV             BQWC119

      PPSPPH=YL2/Y*PYPPH-(VR*PYRPP+VTH*PYTPP+VPH*PYPPP)*YLV             BQWC120

      PNPX=-(2.*U*UX2-YT2*(U-2.*X)+(YT4*(U-2.*X)+4.*YL2*UX*UX2)/RAD)/D2 BQWC121

      PNPY=2.*X*UX*(-YT2+(YT4+2.*YL2*UX2)/RAD)/(D2*Y)                   BQWC122

      PNPZ=I*X*(-2.*UX2-YT2+YT4/RAD)/D2                                 BQWC123

      NNP=N2-(2.*X*PNPX+Y*PNPY+Z*PNPZ)                                  BQWC124

      PNPR =PNPX*PXPR +PNPY*PYPR +PNPZ*PZPR +PNPPS*PPSPR                BQWC125

      PNPTH=PNPX*PXPTH+PNPY*PYPTH+PNPZ*PZPTH+PNPPS*PPSPTH               BQWC126

      PNPPH=PNPX*PXPPH+PNPY*PYPPH+PNPZ*PZPPH+PNPPS*PPSPPH               BQWC127

      PNPVR =PNPPS*(VR *YL2-VDOTY*YR )/V2                               BQWC128

      PNPVTH=PNPPS*(VTH*YL2-VDOTY*YTH)/V2                               BQWC129

      PNPVPH=PNPPS*(VPH*YL2-VDOTY*YPH)/V2                               BQWC130

      PNPT=PNPX*PXPT                                                    BQWC131

      SPACE=REAL(N2).EQ.1..AND.ABS(AIMAG(N2)).LT.ABSLIM                 BQWC132

      POLAR=-I*SQRT(V2)*(-YT2+RAD)/(2.*VDOTY*UX)                        BQWC133

      GAM=(-YT2+RAD)/(2.*UX)                                            BQWC134

      LPOLAR=I*X*SQRT(YT2)/(UX*(U+GAM))                                 BQWC135

      KAY2=OM2/C2*N2                                                    BQWC136

      IF(RSTART.EQ.0.) GO TO 3                                          BQWC137

      SCALE=SQRT(REAL(KAY2)/K2)                                         BQWC138

      KR =SCALE*KR                                                      BQWC139

      KTH=SCALE*KTH                                                     BQWC140

      KPH=SCALE*KPH                                                     BQWC141

  3   CONTINUE                                                          BQWC142

      H=.5*(C2*K2/OM2-N2)                                               BQWC143

      PHPT =-PNPT                                                       BQWC144

      PHPR =-PNPR                                                       BQWC145

      PHPTH=-PNPTH                                                      BQWC146

      PHPPH=-PNPPH                                                      BQWC147

      PHPOM=-NNP/OM                                                     BQWC148

      PHPKR =C2/OM2*KR -C/OM*PNPVR                                      BQWC149

      PHPKTH=C2/OM2*KTH-C/OM*PNPVTH                                     BQWC150

      PHPKPH=C2/OM2*KPH-C/OM*PNPVPH                                     BQWC151

      KPHPK=N2                                                          BQWC152

      RETURN                                                            BQWC153

         END                                                            BQWC154

      SUBROUTINE BQWFNC                                                 BQNC001

C********* CALCULATES A HAMILTONIAN H                                   BQNC002

C********* (= BOOKER QUARTIC FOR VERTICAL INCIDENCE, S=0, C=1)          BQNC003

C********* AND ITS PARTIAL DERIVATIVES WITH RESPECT TO                  BQNC004

C********* TIME, R, THETA, PHI, OMEGA, KR, KTHETA, AND KPHI.            BQNC005

C********* WITH FIELD, NO COLLISIONS                                    BQNC006

      COMMON /CONST/ PI,PIT2,PID2,DEGS,RADIAN,K,C,LOGTEN                BQNC007

      COMMON /RIN/ MODRIN(3),COLL,FIELD,SPACE,KAY2,KAY2I,               BQNC008

     1             H,HI,PHPT,PHPTI,PHPR,PHPRI,PHPTH,PHPTHI,PHPPH,PHPPHI,BQNC009

     2             PHPOM,PHPOMI,PHPKR,PHPKRI,PHPKTH,PHPKTI,PHPKPH,PHPKPIBQNC010

     3            ,KPHPK,KPHPKI,POLAR,POLARI,LPOLAR,LPOLRI,SGN          BQNC011

      COMMON /XX/ MODX(2),X,PXPR,PXPTH,PXPPH,PXPT,HMAX                  BQNC012

      COMMON /YY/ MODY,Y,PYPR,PYPTH,PYPPH,YR,PYRPR,PYRPT,PYRPP,YTH,PYTPRBQNC013

     1           ,PYTPT,PYTPP,YPH,PYPPR,PYPPT,PYPPP                     BQNC014

      COMMON /ZZ/ MODZ,Z(4)                                             BQNC015

      COMMON /RK/ N,STEP,MODE,E1MAX,E1MIN,E2MAX,E2MIN,FACT,RSTART       BQNC016

      COMMON R,TH,PH,KR,KTH,KPH   /WW/ ID(10),W0,W(400)                 BQNC017

      EQUIVALENCE (RAY,W(1)),(F,W(6))                                   BQNC018

      LOGICAL SPACE                                                     BQNC019

      REAL N2,NNP,LPOLAR,LPOLRI,KR,KTH,KPH,K2,KDOTY,K4,KDOTY2,          BQNC020 

     1     KPHPK,KPHPKI,KAY2,KAY2I                                      BQNC021

      DATA (MODRIN=8HBOOKER Q,8HUARTIC, ,8HS=0, C=1),(COLL=0.),         BQNC022

     1     (FIELD=1.),(KAY2I=0.),(HI=0.),(PHPTI=0.),(PHPRI=0.),         BQNC023

     2     (PHPTHI=0.),(PHPPHI=0.),(PHPOMI=0.),(PHPKRI=0.),(PHPKTI=0.), BQNC024

     3     (PHPKPI=0.),(KPHPKI=0.),(POLAR=0.),(LPOLAR=0.),              BQNC025

     4     (X=0.),(PXPR=0.),(PXPTH=0.),(PXPPH=0.),(PXPT=0.),            BQNC026

     5     (Y=0.),(PYPR=0.),(PYPTH=0.),(PYPPH=0.),(YR=0.),(PYRPR=0.),   BQNC027

     6     (PYRPT=0.),(PYRPP=0.),(YTH=0.),(PYTPR=0.),(PYTPT=0.),        BQNC028

     7     (PYTPP=0.),(YPH=0.),(PYPPR=0.),(PYPPT=0.),(PYPPP=0.),        BQNC029

     8     (MODZ=1H ),(U=1.),(U2=1.)                                    BQNC030

      ENTRY RINDEX                                                      BQNC031

      OM=PIT2*1.E6*F                                                    BQNC032

      C2=C*C                                                            BQNC033

      K2=KR*KR+KTH*KTH+KPH*KPH                                          BQNC034

      OM2=OM*OM                                                         BQNC035

      CALL ELECTX                                                       BQNC036

      IF(X.LT..1) GO TO 2                                               BQNC037

      K4=K2*K2                                                          BQNC038

      OM4=OM2*OM2                                                       BQNC039

      C4=C2*C2                                                          BQNC040

      CALL MAGY                                                         BQNC041

      Y2=Y*Y                                                            BQNC042

      KDOTY=KR*YR+KTH*YTH+KPH*YPH                                       BQNC043

      KDOTY2=KDOTY*KDOTY                                                BQNC044

      UX=U-X                                                            BQNC045

      UX2=UX*UX                                                         BQNC046

      A=UX*U2-U*Y2                                                      BQNC047

      B=-2.*U*UX2+Y2*(2.*U-X)                                           BQNC048

      ALPHA=A*C4*K4+X*KDOTY2*C4*K2                                      BQNC049

      BETA=B*C2*OM2*K2-X*KDOTY2*C2*OM2                                  BQNC050

      GAMMA=(UX2-Y2)*UX*OM4                                             BQNC051

      H=ALPHA+BETA+GAMMA                                                BQNC052

      PHPX=-U2*C4*K4+KDOTY2*C4*K2+(4.*U*UX-Y2)*C2*OM2*K2-KDOTY2*C2*OM2+ BQNC053

     1 (-3.*UX2+Y2)*OM4                                                 BQNC054

      PHPY2=-U*C4*K4+(2.*U-X)*C2*OM2*K2-UX*OM4                          BQNC055

      PHPKY2   =X*C2*(C2*K2-OM2)                                        BQNC056

      PHPK2=2.*A*C4*K2+X*KDOTY2*C4+B*C2*OM2                             BQNC057

      PHPT=PHPX*PXPT                                                    BQNC058

      PHPR =PHPX*PXPR +PHPY2*2.*Y*PYPR +PHPKY2  *2.*KDOTY*              BQNC059

     1 (KR*PYRPR+KTH*PYTPR+KPH*PYPPR)                                   BQNC060

      PHPTH=PHPX*PXPTH+PHPY2*2.*Y*PYPTH+PHPKY2  *2.*KDOTY*              BQNC061

     1 (KR*PYRPT+KTH*PYTPT+KPH*PYPPT)                                   BQNC062

      PHPPH=PHPX*PXPPH+PHPY2*2.*Y*PYPPH+PHPKY2  *2.*KDOTY*              BQNC063

     1 (KR*PYRPP+KTH*PYTPP+KPH*PYPPP)                                   BQNC064

      PHPOM=(2.*BETA+4.*GAMMA)/OM                                       BQNC065

     1 -2.*PHPX*X/OM-2.*PHPY2*Y2/OM-2.*PHPKY2  *KDOTY2/OM               BQNC066

      PHPKR= 2.*PHPK2*KR +2.*KDOTY*PHPKY2  *YR                          BQNC067

      PHPKTH=2.*PHPK2*KTH+2.*KDOTY*PHPKY2  *YTH                         BQNC068

      PHPKPH=2.*PHPK2*KPH+2.*KDOTY*PHPKY2  *YPH                         BQNC069

      KAY2 = K2   *(-BETA+RAY*SQRT(BETA**2-4.*ALPHA*GAMMA))/(2.*ALPHA)  BQNC070

      IF(RSTART.EQ.0.) GO TO 1                                          BQNC071

      SCALE=SQRT(KAY2/K2)                                               BQNC072

      KR =SCALE*KR                                                      BQNC073

      KTH=SCALE*KTH                                                     BQNC074

      KPH=SCALE*KPH                                                     BQNC075

  1   CONTINUE                                                          BQNC076

      KPHPK=4.*ALPHA+2.*BETA                                            BQNC077

      SPACE=KAY2.EQ.OM2/C2                                              BQNC078

      POLARI=SQRT(K2)*(U+X*OM2/(C2*KAY2-OM2))/KDOTY                     BQNC079

      LPOLRI= SQRT(Y2-KDOTY2/K2)/UX*(1.-C2*KAY2/OM2)                    BQNC080

      RETURN                                                            BQNC081

C          CALCULATES THE REFRACTIVE INDEX AND ITS GRADIENT USING THE   BQNC082

C          APPLETON-HARTREE FORMULA WITH FIELD, NO COLLISIONS           BQNC083

  2   CONTINUE                                                          BQNC084

      VR =C/OM*KR                                                       BQNC085

      VTH=C/OM*KTH                                                      BQNC086

      VPH=C/OM*KPH                                                      BQNC087

      CALL MAGY                                                         BQNC088

      V2=VR**2+VTH**2+VPH**2                                            BQNC089

      VDOTY=VR*YR+VTH*YTH+VPH*YPH                                       BQNC090

      YLV=VDOTY/V2                                                      BQNC091

      YL2=VDOTY**2/V2                                                   BQNC092

      YT2=Y**2-YL2                                                      BQNC093

      YT4=YT2*YT2                                                       BQNC094

      UX=U-X                                                            BQNC095

      UX2=UX*UX                                                         BQNC096

      RAD=RAY*SQRT(YT4+4.*YL2*UX2)                                      BQNC097

      D=2.*UX-YT2+RAD                                                   BQNC098

      D2=D*D                                                            BQNC099

      N2=1.-2.*X*UX/D                                                   BQNC100

      PNPPS=2.*X*UX*(-1.+(YT2-2.*UX2)/RAD)/D2                           BQNC101

      PPSPR= YL2/Y*PYPR -(VR*PYRPR+VTH*PYTPR+VPH*PYPPR)*YLV             BQNC102

      PPSPTH=YL2/Y*PYPTH-(VR*PYRPT+VTH*PYTPT+VPH*PYPPT)*YLV             BQNC103

      PPSPPH=YL2/Y*PYPPH-(VR*PYRPP+VTH*PYTPP+VPH*PYPPP)*YLV             BQNC104

      PNPX=-(2.*UX2-YT2*(U-2.*X)+(YT4*(U-2.*X)+4.*YL2*UX*UX2)/RAD)/D2   BQNC105

      PNPY=2.*X*UX*(-YT2+(YT4+2.*YL2*UX2)/RAD)/(D2*Y)                   BQNC106

      NNP=N2-(2.*X*PNPX+Y*PNPY)                                         BQNC107

      PNPR =PNPX*PXPR +PNPY*PYPR +PNPPS*PPSPR                           BQNC108

      PNPTH=PNPX*PXPTH+PNPY*PYPTH+PNPPS*PPSPTH                          BQNC109

      PNPPH=PNPX*PXPPH+PNPY*PYPPH+PNPPS*PPSPPH                          BQNC110

      PNPVR =PNPPS*(VR *YL2-VDOTY*YR )/V2                               BQNC111

      PNPVTH=PNPPS*(VTH*YL2-VDOTY*YTH)/V2                               BQNC112

      PNPVPH=PNPPS*(VPH*YL2-VDOTY*YPH)/V2                               BQNC113

      PNPT=PNPX*PXPT                                                    BQNC114

      SPACE=N2.EQ.1.                                                    BQNC115

      POLARI=SQRT(V2)*(YT2-RAD)/(2.*VDOTY*UX)                           BQNC116

      GAM=(-YT2+RAD)/(2.*UX)                                            BQNC117

      LPOLRI=X*SQRT(YT2)/(UX*(U+GAM))                                   BQNC118

      KAY2=OM2/C2*N2                                                    BQNC119

      IF(RSTART.EQ.0.) GO TO 3                                          BQNC120

      SCALE=SQRT(KAY2/K2)                                               BQNC121

      KR =SCALE*KR                                                      BQNC122

      KTH=SCALE*KTH                                                     BQNC123

      KPH=SCALE*KPH                                                     BQNC124

  3   CONTINUE                                                          BQNC125

      H=.5*(C2*K2/OM2-N2)                                               BQNC126

      PHPT =-PNPT                                                       BQNC127

      PHPR =-PNPR                                                       BQNC128

      PHPTH=-PNPTH                                                      BQNC129

      PHPPH=-PNPPH                                                      BQNC130

      PHPOM=-NNP/OM                                                     BQNC131

      PHPKR =C2/OM2*KR -C/OM*PNPVR                                      BQNC132

      PHPKTH=C2/OM2*KTH-C/OM*PNPVTH                                     BQNC133

      PHPKPH=C2/OM2*KPH-C/OM*PNPVPH                                     BQNC134

      KPHPK=N2                                                          BQNC135

      RETURN                                                            BQNC136

         END                                                            BQNC137-

      SUBROUTINE SWWF                                                   SWWF001

C          CALCULATES THE REFRACTIVE INDEX AND ITS GRADIENT USING THE   SWWF002

C          SEN-WYLLER FORMULA -- WITH FIELD                             SWWF003

C          NEEDS SUBROUTINE FSW AND FUNCTIONS C AND S.                  SWWF004

      COMMON /CONST/ PI,PIT2,PID2,DEGS,RADIAN,K,SEA,LOGTEN              SWWF005

      COMMON /RIN/ MODRIN(3),COLL,FIELD,SPACE,KAY2,H,PHPT,PHPR,PHPTH,   SWWF006

     1             PHPPH,PHPOM,PHPKR,PHPKTH,PHPKPH,KPHPK,POLAR,LPOLAR,  SWWF007

     2             SGN                                                  SWWF008

      COMMON /XX/ MODX(2),X,PXPR,PXPTH,PXPPH,PXPT,HMAX                  SWWF009

      COMMON /YY/ MODY,Y,PYPR,PYPTH,PYPPH,YR,PYRPR,PYRPT,PYRPP,YTH,PYTPRSWWF010

     1           ,PYTPT,PYTPP,YPH,PYPPR,PYPPT,PYPPP                     SWWF011

      COMMON /ZZ/ MODZ,Z,PZPR,PZPTH,PZPPH                               SWWF012

      COMMON /RK/ N,STEP,MODE,E1MAX,E1MIN,E2MAX,E2MIN,FACT,RSTART       SWWF013

      COMMON R,TH,PH,KR,KTH,KPH   /WW/ ID(10),W0,W(400)                 SWWF014

      EQUIVALENCE (RAY,W(1)),(F,W(6))                                   SWWF015

      LOGICAL SPACE                                                     SWWF016

      REAL KR,KTH,KPH,K2                                                SWWF017

      COMPLEX KAY2,H,PHPT,PHPR,PHPTH,PHPPH,PHPOM,PHPKR,PHPKTH,PHPKPH,   SWWF018

     1        KPHPK,POLAR,LPOLAR,I,U,RAD,D,PNPPS,PNPX,PNPY,PNPZ,UX,UX2, SWWF019

     2        ALPHA,BETA,GAMMA,A,B,C,TEMP1,TEMP2,TEMP3,ALPOAL,BEPOBE,   SWWF020

     3        GAPOGA,CB2,N2M1,U2,D2GA,DAL,DBET,DGAM,DADY,DADZ,DBDY,DBDZ,SWWF021

     4        DCDY,DCDZ,DUDZ,DT1DX,DT1DY,DT1DZ,DT1DPS,DT2DX,DT2DY,DT2DZ,SWWF022

     5        DT2DPS,DRADDX,DRADDY,DRADDZ,DRDDPS,DDDX,DDDY,DDDZ,DDDPX,  SWWF023

     6        UPX,N2,PNPR,PNPTH,PNPPH,PNPVR,PNPVTH,PNPVPH,NNP,PNPT      SWWF024

      DATA (MODRIN=8H      SE,8HN-WYLLER,8H FORMULA),(COLL=1.),         SWWF025

     1     (FIELD=1.),(LPOLAR=(0.,0.)),                                 SWWF026

     2     (X=0.),(PXPR=0.),(PXPTH=0.),(PXPPH=0.),(PXPT=0.),            SWWF027

     3     (Y=0.),(PYPR=0.),(PYPTH=0.),(PYPPH=0.),(YR=0.),(PYRPR=0.),   SWWF028

     4     (PYRPT=0.),(PYRPP=0.),(YTH=0.),(PYTPR=0.),(PYTPT=0.),        SWWF029

     5     (PYTPP=0.),(YPH=0.),(PYPPR=0.),(PYPPT=0.),(PYPPP=0.),        SWWF030

     6     (Z=0.),(PZPR=0.),(PZPTH=0.),(PZPPH=0.),                      SWWF031

     7     (I=(0.,1.)),(ABSLIM=1.E-5)                                   SWWF032

      ENTRY RINDEX                                                      SWWF033

      OM=PIT2*1.E6*F                                                    SWWF034

      C2=SEA*SEA                                                        SWWF035

      K2=KR*KR+KTH*KTH+KPH*KPH                                          SWWF036

      OM2=OM*OM                                                         SWWF037

      VR =SEA/OM*KR                                                     SWWF038

      VTH=SEA/OM*KTH                                                    SWWF039

      VPH=SEA/OM*KPH                                                    SWWF040

      CALL ELECTX                                                       SWWF041

      CALL MAGY                                                         SWWF042

      OPY=1.+Y                                                          SWWF043

      OMY=1.-Y                                                          SWWF044

      CALL COLFRZ                                                       SWWF045

      Z2=Z*Z                                                            SWWF046

      CALL FSW(1./Z,ALPHA,DAL)                                          SWWF047

      ALPOAL=DAL/ALPHA                                                  SWWF048

      CALL FSW(OMY/Z,BETA,DBET)                                         SWWF049

      BEPOBE=DBET/BETA                                                  SWWF050

      CALL FSW(OPY/Z,GAMMA,DGAM)                                        SWWF051

      GAPOGA=DGAM/GAMMA                                                 SWWF052

      U=Z/ALPHA                                                         SWWF053

      DUDZ=(1.+ALPOAL/Z)/ALPHA                                          SWWF054

      U2=U*U                                                            SWWF055

      UX=U-X                                                            SWWF056

      UPX=U+X                                                           SWWF057

      B=ALPHA/BETA                                                      SWWF058

      DBDY=B*BEPOBE/Z                                                   SWWF059

      DBDZ=-B*(ALPOAL-OMY*BEPOBE)/Z2                                    SWWF060

      C=ALPHA/GAMMA                                                     SWWF061

      DCDY=-C*GAPOGA/Z                                                  SWWF062

      DCDZ=-C*(ALPOAL-OPY*GAPOGA)/Z2                                    SWWF063

      A=.5*(B+C)-1.                                                     SWWF064

      DADY=.5*(DBDY+DCDY)                                               SWWF065

      DADZ=.5*(DBDZ+DCDZ)                                               SWWF066

      TEMP3=(1.-B*C)*U2+A*U*UPX                                         SWWF067

      V2=     VR**2+VTH**2+VPH**2                                       SWWF068

      VDOTY=VR*YR+VTH*YTH+VPH*YPH                                       SWWF069

      YL2=VDOTY**2/V2                                                   SWWF070

      YT2=Y**2-YL2                                                      SWWF071

      Y2=Y*Y                                                            SWWF072

      S2PSI=YT2/Y2                                                      SWWF073

      C2PSI=YL2/Y2                                                      SWWF074

      UX2=UX*UX                                                         SWWF075

      CB2=(C-B)**2                                                      SWWF076

      TEMP1=TEMP3*S2PSI                                                 SWWF077

      DT1DX= A*U*S2PSI                                                  SWWF078

      DT1DY=(U*UPX*DADY-U2*(B*DCDY+C*DBDY))*S2PSI                       SWWF079

      DT1DZ=(2.*U*DUDZ*(1.-B*C+A)+A*X*DUDZ-U2*(B*DCDZ+C*DBDZ)+U*UPX*DADZSWWF080

     1 )*S2PSI                                                          SWWF081

C             (1/YLYT) D/DPSI(TEMP1)                                    SWWF082

      DT1DPS=2.*TEMP1/YT2                                               SWWF083

      TEMP2=U2*CB2*UX2*C2PSI                                            SWWF084

      DT2DX=-2.*UX*U2*CB2*C2PSI                                         SWWF085

      DT2DY=2.*U2*UX2*C2PSI*(C-B)*(DCDY-DBDY)                           SWWF086

      DT2DZ=2.*U2*UX2*C2PSI*(C-B)*(DCDZ-DBDZ)+2.*TEMP2*(1./U+1./UX)*DUDZSWWF087

C             (1/YLYT) D/DPSI(TEMP2)                                    SWWF088

      DT2DPS=-2.*TEMP2/YL2                                              SWWF089

      RAD=RAY*CSQRT(TEMP1**2+TEMP2)                                     SWWF090

      DRADDX=(TEMP1*DT1DX +.5   *DT2DX )/RAD                            SWWF091

      DRADDY=(TEMP1*DT1DY +   .5*DT2DY )/RAD                            SWWF092

      DRADDZ=(TEMP1*DT1DZ +   .5*DT2DZ )/RAD                            SWWF093

C             (1/YLYT) D/DPSI(RAD)                                      SWWF094

      DRDDPS=(TEMP1*DT1DPS+   .5*DT2DPS)/RAD                            SWWF095

      D=2.*U*UX*(1.+A)-TEMP1+RAD+2.*A*U*X*S2PSI                         SWWF096

      DDDX=-2.*U-DT1DX+DRADDX+2.*A*U*S2PSI                              SWWF097

      DDDY= 2.*U*UX*DADY-DT1DY+DRADDY+2.*U*S2PSI*DADY                   SWWF098

      DDDZ=2.*(1.+A)*DUDZ*(U+UX)+2.*U*UX*DADZ-DT1DZ+DRADDZ+2.*X*S2PSI*  SWWF099

     1 (A*DUDZ+U*DADZ)                                                  SWWF100

C          (1/YLYT) D/DPSI(D)                                           SWWF101

      DDDPS=-DT1DPS+DRDDPS+2.*A*U*X/Y2                                  SWWF102

      N2M1=-2.*X*(UX+U*A*S2PSI)/D                                       SWWF103

      N2=1.+N2M1                                                        SWWF104

C               N D/DX(N)                                               SWWF105

      PNPX=-(UX+U*A*S2PSI)*(1.-X*DDDX/D)/D+X/D                          SWWF106

C               N D/DY(N)                                               SWWF107

      PNPY=-X*U*S2PSI/D*DADY-.5*N2M1/D*DDDY                             SWWF108

C               N D/DZ(N)                                               SWWF109

      PNPZ=-X*(1.+A*S2PSI)/D*DUDZ-X*U*S2PSI/D*DADZ-.5*N2M1/D*DDDZ       SWWF110

C             (N/YLYT) D/DPSI(N)                                        SWWF111

      PNPPS=-X*U*A/(D*Y2)     -.5*N2M1/D*DDDPS                          SWWF112

      YLV=VDOTY/V2                                                      SWWF113

C          (YLYT) D/DR(PSI)                                             SWWF114

      PPSPR=YL2/Y*PYPR-(VR*PYRPR+VTH*PYTPR+VPH*PYPPR)*YLV               SWWF115

C          (YLYT) D/DTHETA(PSI)                                         SWWF116

      PPSPTH=YL2/Y*PYPTH-(VR*PYRPT+VTH*PYTPT+VPH*PYPPT)*YLV             SWWF117

C          (YLYT) D/DPHI(PSI)                                           SWWF118

      PPSPPH=YL2/Y*PYPPH-(VR*PYRPP+VTH*PYTPP+VPH*PYPPP)*YLV             SWWF119

      PNPR=PNPX*PXPR+PNPY*PYPR+PNPZ*PZPR+PNPPS*PPSPR                    SWWF120

      PNPTH=PNPX*PXPTH+PNPY*PYPTH+PNPZ*PZPTH+PNPPS*PPSPTH               SWWF121

      PNPPH=PNPX*PXPPH+PNPY*PYPPH+PNPZ*PZPPH+PNPPS*PPSPPH               SWWF122

      PNPVR=PNPPS*(VR*YL2/V2-YLV*YR)                                    SWWF123

      PNPVTH=PNPPS*(VTH*YL2/V2-YLV*YTH)                                 SWWF124

      PNPVPH=PNPPS*(VPH*YL2/V2-YLV*YPH)                                 SWWF125

      NNP=N2-(2.*X*PNPX+Y*PNPY+Z*PNPZ)                                  SWWF126

      PNPT=PNPX*PXPT                                                    SWWF127

      POLAR=I*(TEMP1-RAD)*Y*SQRT(V2)/(U*UX*(C-B)*VDOTY)                 SWWF128

      COSPSI=VDOTY/(Y*SQRT(V2))                                         SWWF129

      LPOLAR=(.5*I*(C-B)*POLAR+A*COSPSI)*SQRT(S2PSI)/                   SWWF130

     1 (POLAR*(UX*(1.+.5*I*(C-B)*COSPSI*POLAR)+A*(U-X*C2PSI)))          SWWF131

      SPACE=REAL(N2).EQ.1..AND.ABS(AIMAG(N2)).LT.ABSLIM                 SWWF132

      KAY2=OM2/C2*N2                                                    SWWF133

      IF(RSTART.EQ.0.) GO TO 1                                          SWWF134

      SCALE=SQRT(REAL(KAY2)/K2)                                         SWWF135

      KR =SCALE*KR                                                      SWWF136

      KTH=SCALE*KTH                                                     SWWF137

      KPH=SCALE*KPH                                                     SWWF138

  1   CONTINUE                                                          SWWF139

C********* CALCULATES A HAMILTONIAN H                                   SWWF140

      H=.5*(C2*K2/OM2-N2)                                               SWWF141

C********* AND ITS PARTIAL DERIVATIVES WITH RESPECT TO                  SWWF142

C********* TIME, R, THETA, PHI, OMEGA, KR, KTHETA, AND KPHI.            SWWF143

      PHPT =-PNPT                                                       SWWF144

      PHPR =-PNPR                                                       SWWF145

      PHPTH=-PNPTH                                                      SWWF146

      PHPPH=-PNPPH                                                      SWWF147

      PHPOM=-NNP/OM                                                     SWWF148

      PHPKR =C2/OM2*KR -SEA/OM*PNPVR                                    SWWF149

      PHPKTH=C2/OM2*KTH-SEA/OM*PNPVTH                                   SWWF150

      PHPKPH=C2/OM2*KPH-SEA/OM*PNPVPH                                   SWWF151

      KPHPK=N2                                                          SWWF152

      RETURN                                                            SWWF153

         END                                                            SWWF154-

      SUBROUTINE SWNF                                                   SWNF001

C          CALCULATES THE REFRACTIVE INDEX AND ITS GRADIENT USING THE   SWNF002

C          SEN-WYLLER FORMULA -- NO FIELD                               SWNF003

C          NEEDS SUBROUTINES FGSW AND FSW AND FUNCTIONS C AND S.        SWNF004

      COMMON /CONST/ PI,PIT2,PID2,DEGS,RADIAN,K,C,LOGTEN                SWNF005

      COMMON /RIN/ MODRIN(3),COLL,FIELD,SPACE,KAY2,H,PHPT,PHPR,PHPTH,   SWNF006

     1             PHPPH,PHPOM,PHPKR,PHPKTH,PHPKPH,KPHPK,POLAR,LPOLAR,  SWNF007

     2             SGN                                                  SWNF008

      COMMON /XX/ MODX(2),X,PXPR,PXPTH,PXPPH,PXPT,HMAX                  SWNF009

      COMMON /YY/ MODY,Y(16)                                            SWNF010

      COMMON /ZZ/ MODZ,Z,PZPR,PZPTH,PZPPH                               SWNF011

      COMMON /RK/ N,STEP,MODE,E1MAX,E1MIN,E2MAX,E2MIN,FACT,RSTART       SWNF012

      COMMON R,TH,PH,KR,KTH,KPH   /WW/ ID(10),W0,W(400)                 SWNF013

      EQUIVALENCE (RAY,W(1)),(F,W(6))                                   SWNF014

      LOGICAL SPACE                                                     SWNF015

      REAL KR,KTH,KPH,K2                                                SWNF016

      COMPLEX KAY2,H,PHPT,PHPR,PHPTH,PHPPH,PHPOM,PHPKR,PHPKTH,PHPKPH,   SWNF017

     1        KPHPK,POLAR,LPOLAR,PNPX,PNPZ,F1,DF,G1,DG1,                SWNF018

     2        N2,PNPR,PNPTH,PNPPH,PNPVR,PNPVTH,PNPVPH,NNP,PNPT          SWNF019

      DATA (MODRIN=8H      SE,8HN-WYLLER,8H FORMULA),(COLL=1.),         SWNF020

     1     (FIELD=0.),(POLAR=(0.,1.)),(LPOLAR=(0.,0.)),                 SWNF021

     2     (X=0.),(PXPR=0.),(PXPTH=0.),(PXPPH=0.),(PXPT=0.),            SWNF022

     3     (MODY=1H ),                                                  SWNF023

     4     (Z=0.),(PZPR=0.),(PZPTH=0.),(PZPPH=0.),                      SWNF024

     5     (ABSLIM=1.E-5),(PNPVR=0.),(PNPVTH=0.),(PNPVPH=0.)            SWNF025

      ENTRY RINDEX                                                      SWNF026

      OM=PIT2*1.E6*F                                                    SWNF027

      C2=C*C                                                            SWNF028

      K2=KR*KR+KTH*KTH+KPH*KPH                                          SWNF029

      OM2=OM*OM                                                         SWNF030

      VR =C/OM*KR                                                       SWNF031

      VTH=C/OM*KTH                                                      SWNF032

      VPH=C/OM*KPH                                                      SWNF033

      CALL ELECTX                                                       SWNF034

      CALL COLFRZ                                                       SWNF035

      CALL FGSW(1./Z,F1,DF1,G1,DG1)                                     SWNF036

      N2=1.-X*G1                                                        SWNF037

      PNPX=-.5*G1                                                       SWNF038

      PNPZ=.5*X*DG1/Z**2                                                SWNF039

      PNPR=PNPX*PXPR+PNPZ*PZPR                                          SWNF040

      PNPTH=PNPX*PXPTH+PNPZ*PZPTH                                       SWNF041

      PNPPH=PNPX*PXPPH+PNPZ*PZPPH                                       SWNF042

      NNP=N2-(2.*X*PNPX+Z*PNPZ)                                         SWNF043

      PNPT=PNPX*PXPT                                                    SWNF044

      SPACE=REAL(N2).EQ.1..AND.ABS(AIMAG(N2)).LT.ABSLIM                 SWNF045

      KAY2=OM2/C2*N2                                                    SWNF046

      IF(RSTART.EQ.0.) GO TO 1                                          SWNF047

      SCALE=SQRT(REAL(KAY2)/K2)                                         SWNF048

      KR =SCALE*KR                                                      SWNF049

      KTH=SCALE*KTH                                                     SWNF050

      KPH=SCALE*KPH                                                     SWNF051

  1   CONTINUE                                                          SWNF052

C********* CALCULATES A HAMILTONIAN H                                   SWNF053

      H=.5*(C2*K2/OM2-N2)                                               SWNF054

C********* AND ITS PARTIAL DERIVATIVES WITH RESPECT TO                  SWNF055

C********* TIME, R, THETA, PHI, OMEGA, KR, KTHETA, AND KPHI.            SWNF056

      PHPT =-PNPT                                                       SWNF057

      PHPR =-PNPR                                                       SWNF058

      PHPTH=-PNPTH                                                      SWNF059

      PHPPH=-PNPPH                                                      SWNF060

      PHPOM=-NNP/OM                                                     SWNF061

      PHPKR =C2/OM2*KR                                                  SWNF062

      PHPKTH=C2/OM2*KTH                                                 SWNF063

      PHPKPH=C2/OM2*KPH                                                 SWNF064

      KPHPK=N2                                                          SWNF065

      RETURN                                                            SWNF066

         END                                                            SWNF067-

      SUBROUTINE FGSW (X,F,DF,G,DG)                                     FGSW001

      COMPLEX F,DF,G,DG                                                 FGSW002

      CALL FSW (X,F,DF)                                                 FGSW003

      IF(ABS(X).GT.50.) GO TO 1                                         FGSW004

      G=X*F                                                             FGSW005

      DG=F+X*DF                                                         FGSW006

      RETURN                                                            FGSW007

    1 X2=X*X                                                            FGSW008

      X3=X2*X                                                           FGSW009

      T2=2.*X2                                                          FGSW010

      T3=3.*X2                                                          FGSW011

      T4=4.*X2                                                          FGSW012

      T8=8.*X2                                                          FGSW013

      T12=12.*X2                                                        FGSW014

      T16=16.*X2                                                        FGSW015

      G=CMPLX(1.-35./T4*(1.-99./T4*(1.-195./T4*(1.-323./T4)))/T4,       FGSW016

     12.5*(1.-63./T4*(1.-143./T4*(1.-255./T4*(1.-399./T4))))/X)         FGSW017

      DG=.5*CMPLX(35.*(1.-99./T2*(1.-585./T8*(1.-323./T3*(1.-2415./T16))FGSW018

     1 ))/X3,                                                           FGSW019

     2-5.*(1.-189./T4*(1.-715./T12*(1.-357./T4*(1.-513./T4))))/X2)      FGSW020

      RETURN                                                            FGSW021

         END                                                            FGSW 22-

      SUBROUTINE FSW (Z,F,DF)                                           FSW 001

C        F(Z) = Z*C3/2(Z) + 2.5*I*C5/2(Z)  AND  DF(Z) = DF/DZ           FSW 002

C        WHERE THE INPUT Z IS REAL AND THE OUTPUT F AND DF ARE COMPLEX. FSW 003

C       NEEDS THE SUBPROGRAMS FOR THE FRESNEL INTEGRAL FUNCTIONS S AND CFSW 004

      DIMENSION A(10),B(10),D(10)                                       FSW 005

      COMPLEX F,DF,C1,C2,C3,C8,W,TEMP,I                                 FSW 006

      DATA (I=(0.,1.)),(PI=3.1415926536),(A3=1.333333333)               FSW 007

      DATA (C2=(1.,1.)),(C3=(1.,-1.)),(C4=.79788456  ),(C6=1.333333333) FSW 008

C        C4=SQRT(2./PI)                                                 FSW 009

      DATA(A=.36230845E-02,.29579186E+00,.23193588E+01,.91355870E+01,   FSW 010

     1.25856287E+02,.60488560E+02,.12562218E+03,.24214980E+03,          FSW 011

     2.44918106E+03,.84244774E+03),                                     FSW 012

     3(B=.16747479E-02,.84796280E-01,.25285001E+00,.22665867E+00,       FSW 013

     4.83871933E-01,.13811875E-01,.98017417E-03,.26299148E-04,          FSW 014

     5.19761006E-06,.18781476E-09),                                     FSW 015

     6(D=.10080653E-03,.46117941E-01,.38507643E+00,.68507885E+00,       FSW 016

     7.42648105E+00,.10742102E+00,.10985920E-01,.40924533E-03,          FSW 017

     8.41881263E-05,.54513142E-08),(G=1.5045055)                        FSW 018

      C1=2./3.*I                                                        FSW 019

      C8=C2*A3*SQRT(PI/2.)                                              FSW 020

      X=Z                                                               FSW 021

      X2=X*X                                                            FSW 022

      X3=X2*X                                                           FSW 023

      IF(ABS(X).GT.50.) GO TO 500                                       FSW 024

      IF(ABS(X).GT.6.) GO TO 1                                          FSW 025

      IF(ABS(X).LT..05)GO TO 200                                        FSW 026

C                              FRESNEL                                  FSW 027

      IF(X.GT.0.) GO TO 300                                             FSW 028

100    Y=C4*SQRT(-X)                                                    FSW 029

      X2=X*X                                                            FSW 030

      W=(COS(X)+I*SIN(X))*(1.-C3*(C(Y)+I*S(Y)))                         FSW 031

      F  =C1+C6*(X+C3*X*X/Y*W)                                          FSW 032

      DF=A3*CMPLX(1.,X)+CMPLX(1.5,X)*A3*C3*X/Y*W                        FSW 033

       RETURN                                                           FSW 034

300   Y=C4*SQRT(X)                                                      FSW 035

      X2=X*X                                                            FSW 036

      W=(COS(X)+I*SIN(X))*(1.-C2*(C(Y)-I*S(Y)))                         FSW 037

      F  =C1+C6*(X-C2*X*X/Y*W)                                          FSW 038

      DF=A3*CMPLX(1.,X)-CMPLX(1.5,X)*A3*C2*X/Y*W                        FSW 039

      RETURN                                                            FSW 040

C                          POWER SERIES                                 FSW 041

  200 X=ABS(Z)                                                          FSW 042

      X2=X*X                                                            FSW 043

      X3=X2*X                                                           FSW 044

      X4=X*X3                                                           FSW 045

      X5=X*X4                                                           FSW 046

      TEMP=-C8* SQRT(X)*CEXP(I*X)                                       FSW 047

      F=CMPLX(4./3.*X-16./9.*X3+64./315.*X5,2./3.+8./3.*X2-32./45.*X4)  FSW 048

     1   +TEMP*X                                                        FSW 049

      DF=CMPLX(4./3.-16./3.*X2+64./63.*X4,16./3.*X-128./45.*X3          FSW 050

     1   +256./945.*X5)                                                 FSW 051

     2   +TEMP*CMPLX(1.5,X)                                             FSW 052

      IF(Z.GE.0.) RETURN                                                FSW 053

      F=-CONJG(F)                                                       FSW 054

      DF=CONJG(DF)                                                      FSW 055

      RETURN                                                            FSW 056

C                               HERMITE                                 FSW 057

    1 XQ = X**2                                                         FSW 058

      X2=XQ                                                             FSW 059

      FR = 0.                                                           FSW 060

      FI = 0.                                                           FSW 061

      DFR = 0.                                                          FSW 062

      DFI = 0.                                                          FSW 063

      DO 2 J = 1,10                                                     FSW 064

      SS = A(J) + XQ                                                    FSW 065

      SB = B(J)/SS                                                      FSW 066

      SD = D(J)/SS                                                      FSW 067

      FR = FR + SB                                                      FSW 068

      FI = FI + SD                                                      FSW 069

      DFR = DFR + SB/SS                                                 FSW 070

    2 DFI = DFI + SD/SS                                                 FSW 071

      F = CMPLX(X*FR,FI)*G                                              FSW 072

      DF = G*(FR - 2.*X*CMPLX(X*DFR,DFI))                               FSW 073

      RETURN                                                            FSW 074

C                            ASYMPTOTIC                                 FSW 075

  500 X2=X*X                                                            FSW 076

      X3=X2*X                                                           FSW 077

      X4=X3*X                                                           FSW 078

      X5=X4*X                                                           FSW 079

      T2=2.*X2                                                          FSW 080

      T3=3.*X2                                                          FSW 081

      T4=4.*X2                                                          FSW 082

      T8=8.*X2                                                          FSW 083

      T16=16.*X2                                                        FSW 084

      T28=28.*X2                                                        FSW 085

      F=CMPLX((1.-35./T4*(1.-99./T4*(1.-195./T4*(1.-323./T4))))/X       FSW 086

     1,5.*(1.-63./T4*(1.-143./T4*(1.-255./T4*(1.-399./T4))))/T2)        FSW 087

      DF=-CMPLX((1.-105./T4*(1.-165./T4*(1.-273./T4*(1.-2907./T28))))/X2FSW 088

     1,5.*(1.-63./T2*(1.-429./T8*(1.-255./T3*(1.-1995./T16))))/X3)      FSW 089

      RETURN                                                            FSW 090

         END                                                            FSW  91-

      FUNCTION C(X)                                                     C   001

      DOUBLEPRECISION   PIH, XD, Y, V, A, QZ, QN, Q, Z                  C   002

      DATA (A1=0.3183099),(A2=0.10132),(B1=0.0968),(B2=0.154)           C   003

      PIH = 1.570796326794897                                           C   004

      XA = ABS(X)                                                       C   005

      IF (XA.GT.4.)   GOTO 20                                           C   006

C                                                                       C   007

      XD = X                                                            C   008

      Y = PIH*XD*XD                                                     C   009

      V = Y*Y                                                           C   010

      A = 1.D0                                                          C   011

      Z = A                                                             C   012

      M = 15.*(XA + 1.)                                                 C   013

      DO 10   I = 1, M                                                  C   014

      KZ=2*(I-1)                                                        C   015

      KV=4*(I-1)                                                        C   016

      QZ = KV + 1                                                       C   017

      QN = (KZ + 1)*(KZ + 2)*(KV + 5)                                   C   018

       Q = QZ/QN                                                        C   019

      A = -A*Q*V                                                        C   020

   10 Z = Z + A                                                         C   021

      Z = Z*XD                                                          C   022

      C = Z                                                             C   023

      RETURN                                                            C   024

C                                                                       C   025

   20  W = PIH*X*X                                                      C   026

      XV=XA**4                                                          C   027

      C=0.5+(A1-B1/XV)*SIN(W)/XA-(A2-B2/XV)*COS(W)/XA**3                C   028

      IF (X.LT.0.)   C = -C                                             C   029

      RETURN                                                            C   030

         END                                                            C    31-

      FUNCTION S(X)                                                     S   001

      DOUBLEPRECISION   PIH, XD, Y, V, A, QZ, QN, Q, Z                  S   002

      DATA (A1=0.3183099),(A2=0.10132),(B1=0.0968),(B2=0.154)           S   003

      PIH = 1.570796326794897                                           S   004

C                                                                       S   005

      XA = ABS(X)                                                       S   006

      IF (XA.GT.4.)   GOTO 20                                           S   007

C                                                                       S   008

      XD = X                                                            S   009

      Y = PIH*XD*XD                                                     S   010

      V = Y*Y                                                           S   011

      A = Y/3.D0                                                        S   012

      Z = A                                                             S   013

      M = 15.*(XA + 1.)                                                 S   014

      DO 10   I = 1, M                                                  S   015

      KZ=2*(I-1)                                                        S   016

      KV=4*(I-1)                                                        S   017

      QZ = KV + 3                                                       S   018

      QN = (KZ + 2)*(KZ + 3)*(KV + 7)                                   S   019

       Q = QZ/QN                                                        S   020

       A = -A*Q*V                                                       S   021

   10 Z = Z + A                                                         S   022

      Z = Z*XD                                                          S   023

      S = Z                                                             S   024

      RETURN                                                            S   025

C                                                                       S   026

   20  W = PIH*X*X                                                      S   027

      XV=XA**4                                                          S   028

      S=0.5-(A1-B1/XV)*COS(W)/XA-(A2-B2/XV)*SIN(W)/XA**3                S   029

      IF (X.LT.0.)   S = -S                                             S   030

      RETURN                                                            S   031

         END                                                            S    32-

