      SUBROUTINE AHWFNC                                                 WFNC001

c	modified by leo mcnamara to handle no-field case as well

c-----	last modified 27-mar-87

	implicit real*8 (a-h,o-z)

	character*8 id,ndate,modx,mody,modz,modrin,koll,modsav

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*8 KR,KTH,KPH,K2,KPHPK,KPHPKI,KAY2,KAY2I,N2,NNP,LPOLAR,LPOLRI WFNC017

	real*8 k,logten

c

	data  ipass / 0 /

	data u / 1. /

c

      ENTRY RINDEX                                                      WFNC027

	if(ray.eq.0.d0.and.ipass.eq.0) print*,'  no magnetic field'

	ipass = 1

      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

c	treat no-field case separately

	if(ray.ne.0.) go to 10

	nnp = 1.

	pnpx = -0.5

	polari = 1.

	lpolri = 1.

	n2 = 1. - x

	go to 20

10    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*dSQRT(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

20	continue

      PNPR =PNPX*PXPR

      PNPTH=PNPX*PXPTH

      PNPPH=PNPX*PXPPH

c	bypass for no-field case

	if(ray.eq.0.) go to 30

      PNPR =PNPr      +PNPY*PYPR +PNPPS*PPSPR                           WFNC056

      PNPTH=PNPth     +PNPY*PYPTH+PNPPS*PPSPTH                          WFNC057

      PNPPH=PNPph     +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

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

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

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

30	continue

      PNPT=PNPX*PXPT                                                    WFNC062

      SPACE=N2.EQ.1.                                                    WFNC063

      KAY2=OM2/C2*N2                                                    WFNC067

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

      SCALE=dSQRT(KAY2/K2)                                              WFNC069

      KR =SCALE*KR                                                      WFNC070

      KTH=SCALE*KTH                                                     WFNC071

      KPH=SCALE*KPH                                                     WFNC072

  1   CONTINUE                                                          WFNC073

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

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

	h = k2 / om2

	h = h * c2

	h = 0.5 * ( h - n2 )

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

      PHPKTH=C2/OM2*KTH

      PHPKPH=C2/OM2*KPH

	if(ray.eq.0.) go to 40

      PHPKR =phpKR     -C/OM*PNPVR                                      WFNC083

      PHPKTH=phpKTH    -C/OM*PNPVTH                                     WFNC084

      PHPKPH=phpKPH    -C/OM*PNPVPH                                     WFNC085

40	continue

      KPHPK=N2                                                          WFNC086

      RETURN                                                            WFNC087

         END                                                            WFNC088-

      block data bdahwfnc                                                WFNC001

c	modified by leo mcnamara to handle no-field case as well

c-----	last modified 20-mar-87

	implicit real*8 (a-h,o-z)

	character*8 id,ndate,modx,mody,modz,modrin,koll,modsav

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*8 KR,KTH,KPH,K2,KPHPK,KPHPKI,KAY2,KAY2I,N2,NNP,LPOLAR,LPOLRI WFNC017

	real*8 k,logten

c

	data coll,kay2i,hi,phpti,phpri,phpthi,phpphi,phpomi,phpkri,

     a       phpkti,phpkpi,kphpki,polar,lpolar,x,pxpr,pxpth,pxpph,pxpt,

     a       y,pypr,pypth,pypph,yr,pyrpr,pyrpt,pyrpp,yth,pytpr,pytpt,

     a       pytpp,yph,pyppr,pyppt,pyppp / 35*0.0 /

	data field / 1.0 /

	data modrin,modz / 'appleton','-hartree',' formula',' ' /

c        data modx(1)/ 'IRIpro' /

      data modx(1) / 'DUDNEY' /

       end

