PROGRAM FLAT
C RAY TRACING IN COMPLEX SPACE IN CARTESIAN COORDINATES (WITH SUBROUTINES)
      COMMON H,X1,X2,VH,V1,V2,RAD(2),P,PP,DUMM(7),INDEP 
      COMMON /WW/ ID(10),DUM,W(400) 
      LOGICAL FIRST,PATH
      REAL INDEP
      COMPLEX H,V,P,PP,AZ1,BETA,N2,FACTOR 
     1,AZEND,AZBEG,AZSTEP,ELEND,ELBEG,ELSTEP
     2,PHIE,THETAE,COSANV,COSANH,ERAY,POLAR,COUPLE
      COMPLEX X1,X2,VH,V1,V2
      DIMENSION DIREC(2),TYPE(2),COUPLE(4)
      DATA (DIREC=4H  UP,4HDOWN),(TYPE=1H ,5HEXTRA) 
      EQUIVALENCE (HBEG,W(324)),(HTRANS,W(1)),(HEND,W(325)),(HREC,W(4)) 
     1,(F,W(10)),(FBEG,W(11)),(FEND,W(12)),(FSTEP,W(13))
     2,(DELTA,W(15)),(DLBEG,W(16)),(DLEND,W(17)),(DLSTEP,W(18)) 
     3,(AZ1,W(20)),(AZBEG,W(22)),(AZEND,W(24)),(AZSTEP,W(26)) 
     4,(BETA,W(30)),(ELBEG,W(32)),(ELEND,W(34)),(ELSTEP,W(36))
     5,(N2,W(301)),(FACTOR,W(319)),(DIP,W(152)) 
     6,(MODX,W(249)),(MODY,W(199)),(MODZ,W(149)),(NTYP,W(7))
     7 ,(FIRST,W(300)),(PATH,W(80)),(POLAR,W(317)),(COUPLE,W(379))
      DATA (PI=3.141592654) 
      CALL Q9EXUN 
      NDATE=IDATE(DUM)
      DEGS=180./PI
  10  CALL INPUT
      PHIH=0. 
      THETAH=1.5*PI-DIP 
      COSTHH=COS(THETAH)
      SINTHH=SIN(THETAH)
      HBEG=HTRANS 
      HEND=HREC 
      NFREQ=1 
      IF (FSTEP.NE.0.) NFREQ=(FEND-FBEG)/FSTEP+1.5
      NDELTA=(DLEND-DLBEG)/DLSTEP+1.5 
                        NAZ=(AZEND-AZBEG)/AZSTEP+1.5
                        NBETA=(ELEND-ELBEG)/ELSTEP+1.5
      DO 50 K=1,NFREQ 
      F=FBEG+(K-1)*FSTEP
      DO 45 J=1,NAZ 
      AZ1=AZBEG+(J-1)*AZSTEP
      AZA=AZ1*DEGS
      DO 40 I=1,NBETA 
      BETA=ELBEG+(I-1)*ELSTEP 
      EL=BETA*DEGS
C     DO 35 NTYP=1,2
      DO 30 L=1,NDELTA
      DELTA=DLBEG+(L-1)*DLSTEP
      FACTOR=CEXP(CMPLX(0.,DELTA))
      DEL=DELTA*180./PI 
      H=HBEG
      X1=X2=0.
      P=0.
      PP=0. 
      INDEP=0.
      V1=-CCOS(BETA)*CCOS(AZ1)
      V2= CCOS(BETA)*CSIN(AZ1)
      N2=1. 
      VH=CSQRT(N2-CCOS(BETA)**2) *SIGN(1.,BETA) 
      FIRST=.TRUE.
      CALL RINDEX 
      VH=CSQRT(N2-CCOS(BETA)**2) *SIGN(1.,BETA) 
      PRINT 1,   MODX,MODY,MODZ,NDATE 
  1   FORMAT(1H1,             4(2X,A8)) 
      FCOSPH=F*SIN(BETA)
      PRINT18,ID(1),ID
  18  FORMAT(1X ,A3,2X,R5,9A8)
      PRINT    17, F,AZA
   17 FORMAT (10X,12HFREQUENCY = ,F14.6,38H MHZ, AZIMUTH ANGLE OF TRANSM
     1ISSION = ,C(F14.6,F14.6), 4H DEG) 
      PRINT 19, FCOSPH,EL,DEL 
   19 FORMAT(* F COSPHI = *,F12.4,*, *
     1             ,34HELEVATION ANGLE OF TRANSMISSION = ,
     2C(F14.6,F14.6),14H DEG, DELTA = ,F14.6, 4H DEG/)
C          VERTICAL POLARIZATION
C     PHIE=PI-AZ1 
C     THETAE=P/2.-BETA
      COSANV=COSTHH*VH+SINTHH*V1
C      HORIZONTAL POLARIZATION
C     PHIE=1.5*PI-AZ1 
C     THETAE=PI/2.
      COSANH=-SINTHH*CSIN(AZ1)
      IF(CABS((V1**2+V2**2)/N2).GT.1.E-10)
     1COSANH=-SINTHH*V1/CSQRT(N2-VH**2) 
C          VERTICAL POLARIZATION
      ERAY=(COSANV+POLAR*COSANH)/(1.-POLAR**2)
      COUPLE(1)=ERAY
      AMPV=20.*ALOG10(CABS(ERAY)) 
      PHASEV=180./PI*CANG(ERAY) 
C      HORIZONTAL POLARIZATION
      ERAY=(COSANH-POLAR*COSANV)/(1.-POLAR**2)
      COUPLE(2)=ERAY
      AMPH=20.*ALOG10(CABS(ERAY)) 
      PHASEH=180./PI*CANG(ERAY) 
      PRINT 7,AMPV,PHASEV,AMPH,PHASEH 
  7   FORMAT( 
     1* ADD*,F10.3,* DB AND*,F10.3,* DEGREES FOR VERTICAL POLARIZATION*/
     2* AND*,F10.3,* DB AND*,F10.3,* DEGREES FOR HORIZONTAL POLARIZATION
     3*/) 
      CALL HEDING 
      NDIREC=2.-SIGN(.5,VH) 
      PRINT 9,DIREC(NDIREC),TYPE(N   TYP) 
  9   FORMAT(1X,A4,*GOING *,A5,*ORDINARY AT TRANSMITTER*) 
      CALL PRINT
      CALL TRACE
      NDIREC=2.-SIGN(.5,VH) 
      NRAYTYP=2.-SIGN(.5,RAD(2))
      PRINT 8,DIREC(NDIREC),TYPE(NRAYTYP) 
  8   FORMAT(1X,A4,*GOING *,A5,*ORDINARY AT RECEIVER*)
      NPUNTYP=2*NTYP-SIGN(.5,RAD(2))
      CALL PRINT
C          VERTICAL POLARIZATION
      COSANV=COSTHH*VH+SINTHH*V1
C      HORIZONTAL POLARIZATION
      COSANH=-SINTHH*CSIN(AZ1)
      IF(CABS((V1**2+V2**2)/N2).GT.1.E-10)
     1COSANH=-SINTHH*V1/CSQRT(N2-VH**2) 
C          VERTICAL POLARIZATION
      ERAY=(COSANV+POLAR*COSANH)/(1.-POLAR**2)
      COUPLE(3)=ERAY
      AMPV=20.*ALOG10(CABS(ERAY)) 
      PHASEV=180./PI*CANG(ERAY) 
C      HORIZONTAL POLARIZATION
      ERAY=(COSANH-POLAR*COSANV)/(1.-POLAR**2)
      COUPLE(4)=ERAY
      AMPH=20.*ALOG10(CABS(ERAY)) 
      PHASEH=180./PI*CANG(ERAY) 
      CALL RAYSET 
      PRINT 7,AMPV,PHASEV,AMPH,PHASEH 
      IF(PATH) CALL RAY END 
      IF(PATH) CALL RAD END 
  30  CONTINUE
  35  CONTINUE
   40 CONTINUE
   45 CONTINUE
   50 CONTINUE
      GO TO 10
         END                                                                  - 
      SUBROUTINE INPUT
      COMMON /WW/ ID(10),DUM,W(400) 
      LOGICAL DEG,INTEG, CYCLE,DIST,P,PINT
C     EQUIVALENCE(NW,W),(EARTHR,W(9)) 
      EQUIVALENCE(NW,W),(EARTHR,W(1)) 
      DATA (PI=3.141592654) 
      DIMENSION NW(400),P(400),PINT(400)
      CALL Q9EXUN 
      READ    1, ID 
  1   FORMAT (10A8) 
      IF(EOF,60) 11,2 
  2   READ    3, K,W(K),DEG,DIST,CYCLE,INTEG
  3   FORMAT (I3,E14.7,4L1) 
      IF (K.LT.0.OR.K.GT.299)  GO TO 9
      IF(K.EQ.0) GO TO 4
      PINT(K)=INTEG 
      P(K)=.NOT.INTEG 
      IF (DEG)  W(K)=W(K)*PI/180. 
      IF (INTEG)NW(K)=W(K)
      IF (CYCLE)  W(K)=W(K)*2.0*PI
      IF (DIST)  W(K)=W(K)/EARTHR 
      GO TO 2 
  4   PRINT 5,ID(1),ID
  5   FORMAT(1H1,1X,A3,2X,R5,9A8) 
      DO  8 K=1,299 
      IF(PINT(K)) PRINT 6,K,NW(K) 
  6   FORMAT(1X,I5,24X,I17) 
      IF(P(K)) PRINT 7,K,W(K) 
  7   FORMAT(1X,I5,2X,E17.10,5X,I17,5X,A8,5X,O16) 
      IF (W(K).NE.0..AND..NOT.P(K).AND..NOT.PINT(K)) PRINT 7,K,W(K) 
     1,NW(K),W(K),W(K)
  8   CONTINUE
      CALL R9EXUN 
      RETURN
  9   PRINT 10,K
  10  FORMAT(* W-ARRAY INPUT DATA OUT OF RANGE, K=* ,I5)
      CALL EXIT 
  11  PRINT 12
  12  FORMAT(*1 OUT OF DATA*) 
      CALL EXIT 
         END                                                                  - 
      SUBROUTINE TRACE
      COMMON H    /WW/ ID(10),DUM,W(400)
      LOGICAL PATH,CHOICE,GROUND,SPACE,THERE
      EQUIVALENCE (MAXSTP,W(61)),(NSKIP,W(60)),(PATH,W(80)),(CHOICE,W(49))
     1)),(SPACE,W(327)),(GROUND,W(326)),(HEND,W(325)),(THERE,W(328))
      CALL DERFUN 
      IF(PATH) CALL RAYPLT
      NPRINT=MAXSTP/NSKIP 
      CHOICE=.TRUE. 
      DO 2 J=1,NPRINT 
      DO 1 K=1,NSKIP
      IF (SPACE) CALL REACH 
      IF(PATH) CALL RAYPLT
      IF (THERE ) GO TO 4 
      HOLD=H
      CALL RKAM 
      IF((H-HEND)*(HOLD-HEND).LT.0.) GO TO 3
      IF(PATH) CALL RAYPLT
      IF(PATH) CALL RADPLT
  1   CONTINUE
      PRINT 6 
  6   FORMAT(1X)
  2   CALL PRINT
      NAXSTP=NSKIP*NPRINT 
      PRINT 5, NAXSTP 
  5   FORMAT(26H THIS RAY NEEDS MORE THAN  ,I5,6H STEPS   ) 
      RETURN
  3   CALL BACK UP
   4  CONTINUE
      IF(PATH) CALL RAYPLT
      RETURN
       END                                                                    - 
      SUBROUTINE REACH
      COMMON HR,HI,X1,X2,VH,V1,V2,RAD   ,P,PP,DUMM( 7),T,ST P 
     1   ,DHDT,DX1DT,DX2DT,DVHDT,DV1DT,DV2DT,DRADDT(2),DPDT,DPPDT 
      COMMON /WW/ ID(10),DUM,W(400) 
      EQUIVALENCE           (GROUND,W(326)),(SPACE,W(327)),(CHOICE, 
     1W(49)),(FACTOR,W(319)),(HBEG,W(324)),(HEND,W(325)),(THERE,W(328)) 
      EQUIVALENCE              (PATH,W(80)),(N2,W(301)) 
      COMPLEXH,V,FACTOR,TEMP,P,PP,DHDT,DVDT,DPDT,DPPDT,T2,T3,H0,N2
      COMPLEX X10,X1,X20,X2,TEMP1,DX1DT,TEMP2,DX2DT,VH0,VH,V10,V1,V20,V2
     1 ,V0,DVHDT,DV1DT,DV2DT ,RAD,RAD0
      LOGICAL GROUND,SPACE,CHOICE,CROSS,THERE 
      LOGICAL PATH,FIRST,START
      EQUIVALENCE(H,HR),(REASTP,W(50)),(V,VR),(FIRST,W(300))
     1,(START,W(377)) 
      DATA(CLOSE=.5E-4) 
      IF(HR.EQ.HBEG) GO TO 6
      PRINT 4 
  4   FORMAT(19H LEAVING IONOSPHERE  )
      CALL PRINT
  6   H0=H
      RAD0=RAD
      X10=X1
      X20=X2
      TEMP=DHDT 
      TEMP1=DX1DT 
      TEMP2=DX2DT 
      VH0=VH
      V10=V1
      V20=V2
      V0=CSQRT(V1**2+V2**2+VH**2) 
      T2=DPDT 
      T3=DPPDT
      S1=(HEND-HR)/REAL(TEMP) 
      CROSS=S1.GT.0.
      STEP=REASTP 
      S=STEP
      DO 1 I=1,200
      IF((S-STEP).GT.S1.AND.CROSS) GO TO 3
      H=H0+S*TEMP 
      X1=X10+S*TEMP1
      X2=X20+S*TEMP2
      START=.TRUE.
      CALL DERFUN 
      IF(SPACE) GO TO 1 
      IF(STEP.LT.CLOSE)GO TO 3
      S=S-STEP
      STEP=STEP/10. 
  1   S=S+STEP
  2   FORMAT(14H LOST IN SPACE) 
      PRINT 2 
      CALL PRINT
      CALL EXIT 
  3   IF(CROSS) S=AMIN1(S,S1) 
      H=H0+S*TEMP 
      X1=X10+S*TEMP1
      X2=X20+S*TEMP2
      P=P+T2*S
      PP=PP+T3*S
      RAD=RAD0
      START=.TRUE.
      CALL RINDEX 
      V1=V10
      V2=V20
      VH=CSQRT(N2-V1**2-V2**2)*SIGN(1.,VH0) 
      THERE =S.EQ.S1
      CHOICE=S.NE.0.
      IF(THERE ) RETURN 
      PRINT 5 
  5   FORMAT(1X,19HENTERING IONOSPHERE  ) 
      CALL PRINT
      RETURN
       END                                                                    - 
      SUBROUTINE BACK UP
      COMMON H,DUMM(24),T,ST P,DRDT 
      COMMON /WW/ ID(10),DUM,W(400) 
      EQUIVALENCE (CHOICE,W(49)),(HEND,W(325))
      EQUIVALENCE (STEP1,W(40)),(INTYP,W(41)) 
      LOGICAL CHOICE
      DATA(CLOSE=.5E-4) 
      DO 1 I=1,10 
      IF(DRDT  .EQ.0.) RETURN 
      STEP=(HEND-H)/DRDT
      STEP=SIGN(AMIN1(ABS(ST P  ),ABS(STEP)),STEP)
      IF (ABS(H        -HEND).LT.CLOSE.AND.STEP.LT.1.)  RETURN
      PRINT 2 
  2   FORMAT(1X,6HHOMING  ) 
      CALL PRINT
      NSAVE=INTYP 
      SAVE=STEP1
      INTYP=1 
      STEP1=STEP
      CHOICE=.TRUE. 
      CALL RKAM 
      INTYP=NSAVE 
      STEP1=SAVE
  1   CHOICE=.TRUE. 
      RETURN
       END                                                                    - 
      SUBROUTINE PRINT
      COMMON HR,HI,X1R,X1I,X2R,X2I,VHR,VHI,V1R,V1I,V2R,V2I,RADR,RADI
     1 ,PR,PIM,PPR,PPI,DUMM( 7) 
     2 ,INDEP,STEP,DHDT,DX1DT,DX2DT,DVHDT,DV1DT,DV2DT,DRADDT(2) 
     3 ,DPDT,DPPDT
      COMMON /WW/ ID(10),DUM,W(400) 
      COMPLEX V,N2,ERROR,I,E,P,VSQR 
     1 ,PP,DHDT,DVDT,DPDT,DPPDT,DT,PPRINT 
     2,BETA,EL,AZA,AZ1
      COMPLEX VH,V1,V2,DX1DT,DX2DT,DVHDT,DV1DT,DV2DT,DP 
      EQUIVALENCE (N2,W(301),T3),(T4,W(302))
     1 ,(F,W(10)),(P,PR)
      EQUIVALENCE          (DELTA,W(15)),(BETA,W(30)),(AZ1,W(20)) 
     1,(VH,VHR),(V1,V1R),(V2,V2R) 
     2,(T1,W(373)),(T2,W(374)),(T5,W(375)),(T6,W(376))
     3,(POLARR,W(317)),(POLARI,W(318))
      DATA(ZERO=0.) 
      DATA (PI=3.141592654) 
      DATA (I=0.,1.),(C=.3) 
      DATA (C=.2997925) 
      REAL LN10 
      DATA (LN10=2.3025851) 
      DP=-I*(VH*HI+V1*X1I+V2*X2I) 
      PPRINT=P+DP 
      PPPRINT=PPR+DP*DPPDT/DPDT 
      E=-2.*PI*F*I*PPRINT/C 
      TIME=PPPRINT/C
      AMP=20./LN10*REAL(E)
      PHASE=180./PI*AIMAG(E) +90. 
      PRINT 2,RADR,T5,VHR,V1R,V2R,T3,T1,HR,X1R,X2R,PR,PPR,POLARR
     1,AMP,PHASE,TIME 
     2 ,RADI,T6,VHI,V1I,V2I,T4,T2,HI,X1I,X2I,PIM,PPI,POLARI 
  2   FORMAT(1X,7E6.0,7F10.4,F11.4,F10.4) 
      RETURN
      ENTRY HEDING
      PRINT 1 
  1   FORMAT(3X,3HRAD,1X,5HERROR,4X,2HVH,4X,2HV1,4X,2HV2,4X,2HN2,1X,
     15HERROR, 9X,1HH, 8X,2HX1, 8X,2HX2, 9X,1HP, 8X,2HP,,5X,5HPOLAR 
     2, 7X, 3HAMP, 6X,5HPHASE, 9X,1HT)
      RETURN
       END                                                                    - 
      SUBROUTINE RAYSET 
      COMMON HR,HI,X1R,X1I,X2R,X2I,VH,V1,V2,RAD(2),PR,PIM,PPR,PPI,DM(7) 
     1 ,INDEP,STEP,DHDT,DX1DT,DX2DT,DVHDT,DV1DT,DV2DT,DRADDT(2) 
     2 ,DPDT,DPPDT
      COMMON /WW/ ID(10),DUM,W(400) 
      EQUIVALENCE (AZI,W(20)),(BETA,W(30))
     1 ,(X1R,X1),(X2R,X2) 
      COMPLEX AZI,BETA,EL,AZ
     1,X1,X2
      COMPLEX VH,V1,V2,DX1DT,DX2DT,DVHDT,DV1DT,DV2DT,DP 
      DIMENSION BUFF(20)
      DATA (BUFF=20(0.)),(ZERO=0.)
      DATA (PI=3.141592654) 
      DATA (I=0.,1.),(C=.3) 
      DATA (C=.2997925) 
      REAL LN10 
      DATA (LN10=2.3025851) 
      EQUIVALENCE (F,W(10)),(NTYP,W(7)),(P,PR),(H,HR) 
      COMPLEX I,E,     H,V,P,PP,D 
     1    ,DHDT,DVDT,DPDT,DPPDT,DT,PPRINT,PPUNCH
      ENTRY PUNCH 
      AZ=180./PI*AZI
      EL=180./PI*BETA 
      DP=-I*(VH*HI+V1*X1I+V2*X2I) 
      PPRINT=P+DP 
      PPPRINT=PPR+DP*DPPDT/DPDT 
      TIME=PPPRINT/C
      D=CSQRT(X1**2+X2**2)
      PPUNCH=PPRINT-CCOS(BETA)*(X2R*CCOS(AZI)-X1R*CSIN(AZI))
     1 +D*CCOS(BETA)-D*COS(BETA)
      E=-2.*PI*F*I*PPUNCH/C 
      AMP=20./LN10*REAL(E)
      PHASE=180./PI*AIMAG(E) +90. 
      PHASE=180.-AMOD(180.-PHASE,360.)
      NPUNTYP=2*NTYP-SIGN(.5,RAD(2))
      NGRONK=NPUNTYP*8+5
      DECODE (8,6,ID) T1,T2 
  6   FORMAT(A1,A2,5X)
  4   IF(UNIT,62) 4,5 
  5   CONTINUE
      ENCODE (144,3,BUFF) T1,NGRONK,T2 ,F,W(20),W(21),W(30),W(31),W(15) 
     1,HR,HI,X1R,X1I,X2R,X2I,AMP,PHASE,TIME,(W(J),J=379,386),ZERO 
  3   FORMAT(A1,R1,A2,23A6,A2)
      BUFFER OUT (62,1)(BUFF(1),BUFF(20)) 
      PRINT 1 
  1   FORMAT(100X,*REFERRED TO TRANSMITTER*)
      PRINT 2,AMP,PHASE 
  2   FORMAT(98X,2F11.4)
      RETURN
       END                                                                    - 
      SUBROUTINE RAYPLT 
      COMMON X,Y
      COMMON /WW/ ID(10),DUM,W(400) 
      DIMENSION XA(600),YA(600),IDAR(5) 
     1,TYPE(2)
      DATA(TYPE=3HORD,3HEXT)
      DATA (N=0),(PI=3.141592654),(DID=5HX3798,1H ,1H ,1H ) 
      COMMON /DDID/DID(4),IFL,LU
      COMMON/DD/IN,IOR,IT,IS,IC,ICC,IX,IY 
      COMMON/DDSCALE/XMIN,XMAX,YMIN,YMAX,MINX,MAXX,MINY,MAXY,DUMM(8)
      EQUIVALENCE (F,W(10)),(AZ,W(20)),(BETA,W(30)),(NTYP,W(7)) 
     1 ,(DELTA,W(15)) 
     2,(CPX,W(85)),(CPY,W(86)),(RX1,W(87)),(RY1,W(88)),(RX2,W(89))
     3,(RY2,W(90))
      DATA(MINX=23),(MAXX=1023),(MINY=23),(MAXY=1023) 
      IF(N.EQ.600) RETURN 
      N=N+1 
      XA(N)=X 
      YA(N)=Y 
      RETURN
      ENTRY RAY END 
      IN=1
      DEL=180./PI*DELTA 
      EL=180./PI*BETA 
      AZDEG=180./PI*AZ
      ENCODE (40,1,IDAR) ID(1),F,DEL,AZDEG,EL,TYPE(NTYP)
   1  FORMAT( A8,*,F=*,F5.3,*,D=*,F5.2,*,A=*,F3.0,*,B=*,F3.0,A3)
      XMIN=W(81)$XMAX=W(82)$YMIN=W(83)$YMAX=W(84) 
      CALL DDGRAPH(N,XA,YA,1,1,1,NRJ) 
      ICC=23B 
      CALL DDGRAPH(1,CPX,CPY,1,1,4,NRJ) 
      ICC=46B 
      CALL DDGRAPH(1,RX1,RY1,1,1,4,NRJ) 
      ICC=67B 
      CALL DDGRAPH(1,RX2,RY2,1,1,4,NRJ) 
  3   FORMAT(F 7.3) 
      IX=8$IY=0 
      CALL DDTAB
      ENCODE( 8,3,NUMBER) XMIN
      CALL DDTABNA8(1,NUMBER,1) 
      IX=943$IY=0 
      CALL DDTAB
      ENCODE( 8,3,NUMBER) XMAX
      CALL DDTABNA8(1,NUMBER,1) 
      IX=0$IY= 8$IOR=1
      CALL DDTAB
      ENCODE( 8,3,NUMBER) YMIN
      CALL DDTABNA8(1,NUMBER,1) 
      IX=0$IY=943 
      CALL DDTAB
      ENCODE( 8,3,NUMBER) YMAX
      CALL DDTABNA8(1,NUMBER,1) 
      IX=20$IY=100$IOR=0
      CALL DDTAB
      CALL DDTABNA8(5,IDAR,1) 
      IN=1
      CALL DDBOXDOT(23,1023,23,1023,100)
      CALL DDBOXDOT(23,1023,23,1023,100)
      CALL DDBOXDOT(23,1023,23,1023,100)
      CALL DDBOXDOT(23,1023,23,1023,100)
      CALL DDBOXDOT(23,1023,23,1023,100)
      CALL DDFRAME
      N=0 
      RETURN
         END                                                                  - 
      SUBROUTINE RADPLT 
      COMMON DUMM(6),X,Y,DUMMM(4),RADR,RADI 
      COMMON /WW/ ID(10),DUM,W(400) 
      DIMENSION XA(600),YA(600),IDAR(5),RAD RA(600),RAD IA(600) 
     1,TYPE(2)
      DATA(TYPE=3HORD,3HEXT)
      DATA (N=1),(PI=3.141592654),(DID=5HX3798,1H ,1H ,1H ) 
     1,(XA=0.),(YA=0.),(RADRA=0.),(RADIA=0.)
      COMMON /DDID/DID(4),IFL,LU
      COMMON/DD/IN,IOR,IT,IS,IC,ICC,IX,IY 
      EQUIVALENCE (F,W(10)),(AZ,W(20)),(BETA,W(30)),(NTYP,W(7)) 
     1 ,(DELTA,W(15)) 
      IF(N.EQ.600) RETURN 
      N=N+1 
      XA(N)=X 
      YA(N)=Y 
      RAD RA(N)=RAD R 
      RAD IA(N)=RAD I 
      RETURN
      ENTRY RAD END 
      DEL=180./PI*DELTA 
      EL=180./PI*BETA 
      AZDEG=180./PI*AZ
      ENCODE (40,1,IDAR) ID(1),F,DEL,AZDEG,EL,TYPE(NTYP)
   1  FORMAT( A8,*,F=*,F5.3,*,D=*,F5.2,*,A=*,F3.0,*,B=*,F3.0,A3)
      CALL FILMGRAF(RAD RA,RAD IA,N,IDAR,1,0) 
      CALL FILMGRAF(XA,YA,N,IDAR,1,0) 
      N=1 
      RETURN
         END                                                                  - 
      SUBROUTINE RKAM 
C        NUMERICAL INTEGRATION OF DIFFERENTIAL EQUATIONS
C     COMMON/SHARE/NN,SPACE ,MODE,KKA,E1MAX,E1MIN,E2MAX,E2MIN,FACT          3 
      COMMON Y(25),T,STEP,DYDT(25)
      COMMON /WW/ ID(10),DUM,W(400) 
      EQUIVALENCE(SPACE,W(40)),(MODE,W(41)),(E1MAX,W(42)),(ERATIO,W(43))
     1,(KKA,W(44)),(E2MAX,W(45)),(E2MIN,W(46)),(FACT,W(47)) 
     2,(NN,W(48)),(RSTAR,W(49)) 
      DIMENSION DELY(4,25),BET(25),XV(5),FV(4,25),YU(5,25)
C     TYPE DOUBLE YU                                                        5 
      DOUBLE PRECISION YU 
      LOGICAL RSTAR 
      IF(RSTAR) GO TO 10
      GO TO (1001,2000,2000),MODE                                           6 
C      RUNGE-KUTTA                                                          7 
 1000 LL=1                                                                  8 
 1001 DO 1034 K=1,4 
      DO 1350 I=1,NN
      DELY(K,I)=STEP*FV(MM,I) 
      Z=YU(MM,I)
 1350 Y(I)=Z+BET(K)*DELY(K,I) 
      T=BET(K)*STEP+XV(MM)
      CALL DERFUN 
      DO 1034 I=1,NN
 1034 FV(MM,I)=DYDT(I)
      DO 1039 I=1,NN
      DEL=(DELY(1,I)+2.0*DELY(2,I)+2.0*DELY(3,I)+DELY(4,I))/6.0 
 1039 YU(MM+1,I)=YU(MM,I)+DEL 
      MM=MM+1 
      XV(MM)=XV(MM-1)+STEP
      DO 1400 I=1,NN
 1400 Y(I)=YU(MM,I) 
      T=XV(MM)
      CALL DERFUN 
      GO TO (42,100,100),MODE                                              19 
  100 DO 150 I=1,NN                                                        20 
  150 FV(MM,I)=DYDT(I)
      GO TO (1001,1001,1001,2000),MM
C       ADAMS-MOULTON                                                      22 
 2000 DO 2048 I=1,NN                                                       23 
      DEL=STEP*(55.0*FV(4,I)-59.0*FV(3,I)+37.0*FV(2,I)-9.0*FV(1,I))/24.   24
      Y(I)=YU(4,I)+DEL                                                     25 
 2048 DELY(1,I)=Y(I)
      T=XV(4)+STEP
      CALL DERFUN 
      XV(5)=T 
      DO 2051 I=1,NN
      DEL=STEP*(9.0*DYDT(I)+19.0*FV(4,I)-5.0*FV(3,I)+FV(2,I))/24.0
      YU(5,I)=YU(4,I)+DEL                                                  29 
 2051 Y(I)=YU(5,I)
      CALL DERFUN 
      GO TO (42,42,3000),MODE                                              31 
C        ERROR ANALYSIS                                                    32 
 3000 SSE=0.0 
      DO 3033 I=1,NN
      EPSIL=R*ABS (Y(I)-DELY(1,I))
      GO TO (3301,3307),KKA 
 3301 IF(Y(I)) 3650,3307,3650                                              35 
 3650 EPSIL=EPSIL/ABS (Y(I))
 3307 IF(SSE-EPSIL) 3032,3033,3033                                         37 
 3032 SSE=EPSIL                                                            38 
 3033 CONTINUE
      IF (E1MAX-SSE) 3034,3034,3035 
 3034 IF(ABS (STEP)-E2MIN) 42,42,4340 
 3035 IF(SSE-E1MIN) 3036,42,42                                             41 
 3036 IF(E2MAX-ABS (STEP)) 42,42,5360 
 4340 LL=1                                                                 43 
      MM=1                                                                 44 
      STEP=STEP*FACT
      GO TO 1001                                                           46 
 5360 GO TO (42,5361),LL                                                   47 
 5361 XV(2)=XV(3) 
      XV(3)=XV(5) 
      DO 5363 I=1,NN
      FV(2,I)=FV(3,I) 
      FV(3,I)=DYDT(I) 
      YU(2,I)=YU(3,I)                                                      50 
 5363 YU(3,I)=YU(5,I) 
      STEP=2.0*STEP 
      LL=2
      MM=3
      GO TO 1001
C       EXIT ROUTINE                                                       52 
   42 DO 12 K=1,3 
      XV(K)=XV(K+1) 
      DO 12 I=1,NN
      FV(K,I)=FV(K+1,I) 
   12 YU(K,I)=YU(K+1,I) 
      LL=2
      MM=4
      XV(4)=XV(5) 
      DO 52 I=1,NN
      FV(4,I)=DYDT(I) 
   52 YU(4,I)=YU(5,I) 
      GO TO (70,70,73),MODE 
C 
      ENTRY RSTART
C 
  10  CONTINUE
      RSTAR=.FALSE. 
      ALPHA=T 
      EPM=0.0 
      GO TO (7,9,9),MODE
    7 MM=4
      GO TO 8 
    9 MM=1                                                                 60 
    8 BET(1)=0.5
      BET(2)=0.5
      BET(3)=1.0
      BET(4)=0.0
      STEP=SPACE
      R=19.0/270.0
      XV(MM)=T
      IF(ERATIO.LE.0.) ERATIO=55. 
      E1MIN=E1MAX/ERATIO
      IF (FACT.LE.0.) FACT=0.5
      CALL DERFUN 
      DO 320 I=1,NN 
      FV(MM,I)=DYDT(I)
  320 YU(MM,I)=Y(I) 
      GO TO 1000
   73 E=ABS (XV(4)-ALPHA) 
      IF (E-EPM) 2000,2000,71 
   71 EPM=E                                                                70 
   70 RETURN                                                               71 
      END                                                                     - 
      SUBROUTINE DERFUN 
C        HASELGROVES EQUATION FOR FLAT EARTH
      COMMON H,X1,X2,VH,V1,V2,RAD,P,PP,DUMM(7),INDEP,STEP 
     1 ,DHDT,DX1DT,DX2DT,DVHDT,DV1DT,DV2DT,DRADDT,DPDT,DPPDT
      COMPLEX DRADDT,RAD
      COMMON /WW/ ID(10),DUM,W(400) 
      COMPLEX H,V,P,PP,DHDT,DVDT,DPDT,DPPDT,FACTOR,NDNDH,N2,NNP,TEMP
      COMPLEX X1,X2,VH,V1,V2,DX1DT,DX2DT,DVHDT,DV1DT,DV2DT
      COMPLEX NDNDX1,NDNDX2,NDNDVH,NDNDV1,NDNDV2
      COMPLEX I 
      DATA (I=0.,1.)
      EQUIVALENCE (UX,W(329)),(UX2,W(331)),(YT2,W(333)),(YL2,W(335))
     1 ,(PXPH,W(337)),(PXPX1,W(339)),(PXPX2,W(341)) 
     2 ,(PZPH,W(343)),(PZPX1,W(345)),(PZPX2,W(347)) 
     3 ,(PYLPH,W(349)),(PYLPX1,W(351)),(PYLPX2,W(353))
     4 ,(PYLPVH,W(355)),(PYLPV1,W(357)),(PYLPV2,W(359)) 
     5 ,(PYTPH,W(361)),(PYTPX1,W(363)),(PYTPX2,W(365))
     6 ,(PYTPVH,W(367)),(PYTPV1,W(369)),(PYTPV2,W(371)) 
      COMPLEX DXDT,DYLDT,DYTDT,DZDT,DUDT, 
     1  UX,UX2,YT2,YL2,PYLPH,PYLPX1,PYLPX2,PYLPVH,PYLPV1,PYLPV2,PYTPH 
     2 ,PYTPX1,PYTPX2,PYTPVH,PYTPV1,PYTPV2,PXPH,PXPX1,PXPX2 
     3 ,PZPH,PZPX1,PZPX2
      EQUIVALENCE (FACTOR,W(319)),(NDNDH,W(305)),(N2,W(301)),(NNP,W(303)) 
     1),(SPACE,W(327)),(SPAERR,W(51)) 
      EQUIVALENCE (NDN       DX1,W(307)),(NDNDX2,W(309))
     1 ,(NDNDVH,W(311)),(NDNDV1,W(313)),(NDNDV2,W(315)) 
      LOGICAL SPACE 
      CALL RINDEX 
      DHDT=FACTOR*(VH-NDNDVH) 
      DX1DT=FACTOR*(V1-NDNDV1)
      DX2DT=FACTOR*(V2-NDNDV2)
      DVHDT=FACTOR*NDNDH
      DV1DT=FACTOR*NDNDX1 
      DV2DT=FACTOR*NDNDX2 
      DXDT=PXPH*DHDT+PXPX1*DX1DT+PXPX2*DX2DT
      DYLDT=PYLPH*DHDT+PYLPX1*DX1DT+PYLPX2*DX2DT
     1 +PYLPVH*DVHDT+PYLPV1*DV1DT+PYLPV2*DV2DT
      DYTDT=PYTPH*DHDT+PYTPX1*DX1DT+PYTPX2*DX2DT
     1 +PYTPVH*DVHDT+PYTPV1*DV1DT+PYTPV2*DV2DT
      DZDT=PZPH*DHDT+PZPX1*DX1DT+PZPX2*DX2DT
      DUDT=-I*DZDT
      DRADDT=(2.*YT2*DYTDT+4.*DYLDT*UX2+4.*YL2*UX*(DUDT-DXDT))/RAD
      DPDT=FACTOR*N2
      DPPDT=FACTOR*NNP
      DVDT=CSQRT(DVHDT**2+DV1DT**2+DV2DT**2)
      SPACE=CABS(DVDT).LE.SPAERR
      RETURN
       END                                                                    - 
      SUBROUTINE WFWC 
C                       CALCULATES THE REFRACTIVE INDEX AND ITS GRADIENT
C        INCLUDING THE EARTHS MAGNETIC FIELD AND INCLUDING ELECTRON COLLISIONS. 
C        REWRITTEN OCTOBER 1967 
      COMMON H,TH,PH,VR,VTH,VPH,RADR,RADI 
      EQUIVALENCE (UX,W(329)),(UX2,W(331)),(YT2,W(333)),(YL2,W(335))
     1 ,(PXPR,W(337)),(PZPR,W(343)) 
     2 ,(PYLPR,W(349)),(PYLPTH,W(351)),(PYLPPH,W(353))
     3 ,(PYLPVR,W(355)),(PYLPVT,W(357)),(PYLPVP,W(359)) 
     4 ,(PYTPR,W(361)),(PYTPTH,W(363)),(PYTPPH,W(365))
     5 ,(PYTPVR,W(367)),(PYTPVT,W(369)),(PYTPVP,W(371)) 
      COMPLEX                   PYLPR,PYLPTH,PYLPPH,PYLPVR,PYLPVT,PYLPVP
     1,PYTPR,PYTPTH,PYTPPH,PYTPVR,PYTPVT,PYTPVP,YLV,VDOTY 
      DIMENSION PXPR(3),YR(4),PYRPR(4,3),PZPR(3),A(2) 
      COMMON      /W/ ID(10),DUM,W(400) 
      EQUIVALENCE (NTYP,W(7)),(N2,W(301)),(NNP,W(303)),(PNPR,W(305))
     1 ,(PNPTH,W(307)),(PNPPH,W(309)),(PNPVR,W(311)),(PNPVTH,W(313))
     2 ,(PNPVPH,W(315)),(POLAR,W(317))
      EQUIVALENCE 
     1 (PXPTH,PXPR(2)),(PXPPH,PXPR(3)),(YTH,YR(2)),(YPH,YR(3)),(Y,YR(4))
     2 ,(PYRPTH,PYRPR(5)),(PYRPPH,PYRPR(9)),(PYTHPR,PYRPR(2)),
     3 (PYTHPTH,PYRPR(6)),(PYTHPPH,PYRPR(10)),(PYPHPR,PYRPR(3)),
     4 (PYPHPTH,PYRPR(7)),(PYPHPPH,PYRPR(11)),(PYPR,PYRPR(4)),
     5 (PYPTH,PYRPR(8)),(PYPPH,PYRPR(12)),(PZPTH,PZPR(2)),
     6 (PZPPH,PZPR(3))
      COMPLEX I,U,RAD,D,N2,PNPPS,PNPX,PNPY,PNPZ,POLAR,UX,UX2,D2,RAD2
      COMPLEX NNP,X,PXPR,Z,PZPR,PNPR,YR,PYRPR,PNPTH,PNPPH,PNPVR,PNPVTH
     1 ,PNPVPH,PXPTH,PXPPH,YTH,YPH,H,TH,PH,VR,VTH,VPH,V2,V,COSPS,SINPS
     2 ,YT,YT2,YT4,YL,YL2,PPSPR,PPSPTH,PPSPPH,PYRPTH,PYRPPH,PYTHPR
     3 ,PYTHPTH,PYTHPPH,PYPHPR,PYPHPTH,PYPHPPH,PYPR,PYPTH,PYPPH,PZPTH 
     4 ,PZPPH,VERR,RADERR,CONST 
      DATA (I=0.,1.),(X=2(0.)),(PXPR=6(0.)),(YR=8(0.)),(PYRPR=24(0.)) 
     1 ,(Z=2(0.)),(PZPR=6(0.)),(CONST=.7071067812,.7071067812)
C     CONST=CSQRT(I)
      EQUIVALENCE (START,W(377)),(RAD,RADR),(RAD2,RAD2R,A),(RAD2I,A(2)) 
     1,(PATH,W(80)),(FIRST,W(300)),(VERR,W(373)),(RADERR,W(375))
      LOGICAL PATH,FIRST,START
      ENTRY RINDEX
      CALL ELECTX (X,PXPR)
      CALL MAGY (YR,PYRPR)
      V2=     VR**2+VTH**2+VPH**2 
C     COSPS=(VR*YR+VTH*YTH+VPH*YPH)/(V*Y) 
C     SINPS=CSQRT(1.-COSPS**2)
C     YT=Y*SINPS
C     YT2=YT*YT 
      VDOTY=VR*YR+VTH*YTH+VPH*YPH 
      YLV=VDOTY/V2
      YL2=VDOTY**2/V2 
      YT2=Y**2-YL2
      YT4=YT2*YT2 
C     YL=Y*COSPS
C     YL2=YL*YL 
C 
C           YL*P(YL)/P(R)          ETC. 
      PYLPR=YLV*(VR*PYRPR+VTH*PYTHPR+VPH*PYPHPR)
      PYLPTH=YLV*(VR*PYRPTH+VTH*PYTHPTH+VPH*PYPHPTH)
      PYLPPH=YLV*(VR*PYRPPH+VTH*PYTHPPH+VPH*PYPHPPH)
      PYLPVR=YLV*YR-YL2*VR/V2 
      PYLPVT=YLV*YTH-YL2*VTH/V2 
      PYLPVP=YLV*YPH-YL2*VPH/V2 
C 
C            YT*P(YT)/P(R)     ETC. 
      PYTPR=Y*PYPR-PYLPR
      PYTPTH=Y*PYPTH-PYLPTH 
      PYTPPH=Y*PYPPH-PYLPPH 
      PYTPVR=-PYLPVR
      PYTPVT=-PYLPVT
      PYTPVP=-PYLPVP
      CALL COLFRZ (Z,PZPR)
      U=1.-I*Z
      UX=U-X
      UX2=UX*UX 
      RAD2=YT4+4.*YL2*UX2 
      IF(FIRST) RADI=1.5-NTYP 
      IF(START.OR.FIRST) RAD=CSQRT(-RAD2)*I*SIGN(1.,RADI) 
      FIRST=START=.FALSE. 
      RADERR=RAD**2/RAD2-1. 
      D=2.*U*UX-YT2+RAD 
      D2=D*D
      N2=1.-2.*X*UX/D 
      VERR=V2/N2-1. 
      PNPPS=2.*X*UX   *(-1.+(YT2-2.*UX2)/RAD)/D2
C         WE HAVE DIVIDED PNPPS BY YL SEPT 30, 1969 
C           AND AT THE SAME TIME, MULTIPLIED THE NEXT 3 VARIABLES BY YL 
      PPSPR= YL2/Y*PYPR-(VR*PYRPR+VTH*PYTHPR+VPH*PYPHPR)*YLV
      PPSPTH= YL2/Y*PYPTH-(VR*PYRPTH+VTH*PYTHPTH+VPH*PYPHPTH)*YLV 
      PPSPPH= YL2/Y*PYPPH-(VR*PYRPPH+VTH*PYTHPPH+VPH*PYPHPPH)*YLV 
      PNPX=-(2.*U*UX2-YT2*(U-2.*X)+(YT4*(U-2.*X)+4.*YL2*UX*UX2)/RAD)/D2 
      PNPY=2.*X*UX*(-YT2+(YT4+2.*YL2*UX2)/RAD)/(D2*Y) 
      PNPZ=I*X*(-2.*UX2-YT2+YT4/RAD)/D2 
      PNPR=PNPX*PXPR+PNPY*PYPR+PNPZ*PZPR+PNPPS*PPSPR
      PNPTH=PNPX*PXPTH+PNPY*PYPTH+PNPZ*PZPTH+PNPPS*PPSPTH 
      PNPPH=PNPX*PXPPH+PNPY*PYPPH+PNPZ*PZPPH+PNPPS*PPSPPH 
      PNPVR=PNPPS*(VR*YL2/V2-YLV*YR)
      PNPVTH=PNPPS*(VTH*YL2/V2-YLV*YTH) 
      PNPVPH=PNPPS*(VPH*YL2/V2-YLV*YPH) 
      NNP=N2-(2.*X*PNPX+Y*PNPY+Z*PNPZ)
      V=CSQRT(V2) 
      YL=VDOTY/V
      POLAR=-I*(-YT2+RAD)/(2.*YL*UX)
      RETURN
      END                                                                     - 
      SUBROUTINE NFWC 
C        WITH COLLISIONS, NO FIELD
C            VERSION FOR COMPLEX SPACE
      COMMON      /W/ ID(10),DUM,W(400) 
      COMPLEX U,N2,I,NNP,X,PXPR,Z,PZPR,PNPX,PNPZ,PNPR 
      EQUIVALENCE(N2,W(301)),(NNP,W(303)),(PNPR,W(305)),(FIELD,W(300))
     1 ,(CALLX,W(330)),(CALLZ,W(340)),(X,W(331)),(PXPR,W(333))
     2 ,(Z,W(341)),(PZPR,W(343)),(MODY,W(199))
      DIMENSION PXPR(3),PZPR(3),PNPR(3) 
      DATA(X=2(0.)),(Z=2(0.)),(PXPR=6(0.)),(PZPR=6(0.)),(I=0.,1.) 
      DATA (MODY=8HNO FIELD)
      ENTRY RINDEX
      CALL ELECTX (X,PXPR)
      CALL COLFRZ (Z,PZPR)
      U=1.-I*Z
      N2=1.-X/U 
      PNPX=-1./(2.*U) 
      PNPZ=     -I*X/(2.*U**2)
      NNP=1.+I*X*Z/U**2 
      PNPR=PNPX*PXPR+PNPZ*PZPR
      PNPR(2)=PNPX*PXPR(2)+PNPZ*PZPR(2) 
      PNPR(3)=PNPX*PXPR(3)+PNPZ*PZPR(3) 
      RETURN
      ENTRY MAGY
      RETURN
       END                                                                    -