SUBROUTINE PRINTR(NWHY,CARD) PRIN001
c$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
c note format 1200 changed 17-dec-86 3P2f6.0 --> 2P2f6.0
c$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
c there is an inconsistency in the way POLAR is treated
c This gives 4 values in the print-out e.g PRIN166.
c The group path is the last column printed (no header)
c$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
implicit real*8 (a-h,o-z)
character*8 id,ndate,modx,mody,modz,modrin,koll,modsav
character*8 headr1,headr2,head1,head2,unit,units
character*1 type
C PRINTS OUTPUT AND PUNCHES RAYSETS WHEN REQUESTED PRIN002
DIMENSION G(3,3),G1(3,3),TYPE(3),HEADR1(20),HEADR2(20),UNITS(20), PRIN003
1 HEAD1(20),HEAD2(20),UNIT(20),RPRINT(20),NPR(20) PRIN004
COMMON /CONST/ PI,PIT2,PID2,DEGS,RAD,DUM(3) PRIN005
COMMON /FLG/ NTYP,NEWWR,NEWWP,PENET,LINES,IHOP,HPUNCH PRIN006
COMMON /RIN/ MODRIN(3),COLL,FIELD,SPACE,N2,N2I,PNP(10),POLAR(2), PRIN007
1 LPOLAR(2) PRIN008
COMMON R(20),T /WW/ ID(10),W0,W(400) PRIN009
common / leo / apogee,range
EQUIVALENCE (THETA,R(2)),(PHI,R(3)) PRIN010
EQUIVALENCE (EARTHR,W(2)),(XMTRH,W(3)),(TLAT,W(4)),(TLON,W(5)), PRIN011
1 (F,W(6)),(AZ1,W(10)),(BETA,W(14)),(RCVRH,W(20)),(HOP,W(22)), PRIN012
2 (PLAT,W(24)),(PLON,W(25)),(RAYSET,W(72)) PRIN013
LOGICAL SPACE,NEWWR,NEWWP,PENET PRIN014
REAL*8 N2,N2I PRIN015
COMPLEX*16 PNP,polar,lpolar PRIN016
DATA TYPE / 'X','N','O'/ PRIN017
2,HEADR1(7)/' PHAS'/,HEADR2(7)/'E PATH'/,UNITS(7)/' KM' /, PRIN018
3 HEADR1(8)/' ABSO'/,HEADR2(8)/'RPTION'/,UNITS(8)/' DB' /, PRIN019
4 HEADR1(9)/' DOP'/,HEADR2(9)/'PLER '/,UNITS(9)/' C/S' /, PRIN020
5 HEADR1(10)/' PATH' /,HEADR2(10)/'LENGTH'/,UNITS(10)/' KM' / PRIN023
c
data cos_ell_save / 1.d0/
c magnetic dip angle - r(2) = colatitude. Dipole field
tan_dip = 2.d0 / dtan(r(2))
dipole_dip = datan(tan_dip)
c
CALL RINDEX PRIN024
IF (.NOT.NEWWP) GO TO 10 PRIN025
C********* NEW W ARRAY -- REINITIALIZE PRIN026
NEWWP=.FALSE. PRIN027
SPL=dSIN(PLON-TLON) PRIN028
CPL=dSIN(PID2-(PLON-TLON)) PRIN029
SP=dSIN(PLAT) PRIN030
CP=dSIN(PID2-PLAT) PRIN031
SL=dSIN(TLAT) PRIN032
CL=dSIN(PID2-TLAT) PRIN033
C********* MATRIX TO ROTATE COORDINATES PRIN034
G(1,1)=CPL*SP*CL-CP*SL PRIN035
G(1,2)=SPL*SP PRIN036
G(1,3)=-SL*SP*CPL-CL*CP PRIN037
G(2,1)=-SPL*CL PRIN038
G(2,2)=CPL PRIN039
G(2,3)=SL*SPL PRIN040
G(3,1)=CL*CP*CPL+SP*SL PRIN041
G(3,2)=CP*SPL PRIN042
G(3,3)=-SL*CP*CPL+SP*CL PRIN043
DENM=G(1,1)*G(2,2)*G(3,3)+G(1,2)*G(3,1)*G(2,3)+G(2,1)*G(3,2)*G(1,3PRIN044
1)-G(2,2)*G(3,1)*G(1,3)-G(1,2)*G(2,1)*G(3,3)-G(1,1)*G(3,2)*G(2,3) PRIN045
C********* THE MATRIX G1 IS THE INVERSE OF THE MATRIX G PRIN046
G1(1,1)=(G(2,2)*G(3,3)-G(3,2)*G(2,3))/DENM PRIN047
G1(1,2)=(G(3,2)*G(1,3)-G(1,2)*G(3,3))/DENM PRIN048
G1(1,3)=(G(1,2)*G(2,3)-G(2,2)*G(1,3))/DENM PRIN049
G1(2,1)=(G(3,1)*G(2,3)-G(2,1)*G(3,3))/DENM PRIN050
G1(2,2)=(G(1,1)*G(3,3)-G(3,1)*G(1,3))/DENM PRIN051
G1(2,3)=(G(2,1)*G(1,3)-G(1,1)*G(2,3))/DENM PRIN052
G1(3,1)=(G(2,1)*G(3,2)-G(3,1)*G(2,2))/DENM PRIN053
G1(3,2)=(G(3,1)*G(1,2)-G(1,1)*G(3,2))/DENM PRIN054
G1(3,3)=(G(1,1)*G(2,2)-G(2,1)*G(1,2))/DENM PRIN055
R0=EARTHR+XMTRH PRIN056
C********* CARTESIAN COORDINATES OF TRANSMITTER PRIN057
XR=R0*G(1,1) PRIN058
YR=R0*G(2,1) PRIN059
ZR=R0*G(3,1) PRIN060
CTHR=G(3,1) PRIN061
STHR=dSIN(dACOS (CTHR)) PRIN062
PHIR=dATAN2(YR,XR) PRIN063
ALPH=dATAN2(G(3,2),G(3,3)) PRIN064
C********* PRIN065
NR=6 PRIN066
NP=0 PRIN067
DO 7 NN=7,20 PRIN068
IF (W(NN+50).EQ.0.) GO TO 7 PRIN069
C********* DEPENDENT VARIABLE NUMBER NN IS BEING INTEGRATED PRIN070
C********* NR IS THE NUMBER OF DEPENDENT VARIABLES BEING INTEGRATED PRIN071
NR=NR+1 PRIN072
IF (W(NN+50).NE.2.) GO TO 7 PRIN073
C********* DEPENDENT VARIABLE NUMBER NN IS BEING INTEGRATED AND PRINTED.PRIN074
C********* NP IS THE NUMBER OF DEPENDENT VARIABLES BEING INTEGRATED AND PRIN075
C********* PRINTED PRIN076
NP=NP+1 PRIN077
C********* SAVE THE INDEX OF THE DEPENDENT VARIABLE TO PRINT PRIN078
NPR(NP)=NR PRIN079
HEAD1(NP)=HEADR1(NN) PRIN080
HEAD2(NP)=HEADR2(NN) PRIN081
UNIT(NP)=UNITS(NN) PRIN082
7 CONTINUE PRIN083
NP1=MIN0(NP,3) PRIN084
ccc PDEV=ABSORB=DOPP=0. PRIN085
pdev = 0.d0
absorb = 0.d0
dopp = 0.d0
C********* PRINT COLUMN HEADINGS AT THE BEGINNING OF EACH RAY PRIN086
10 IF (IHOP.NE.0.d0) GO TO 12 PRIN087
write(2,1100)(HEAD1(NN),HEAD2(NN),NN=1,NP1) PRIN088
1100 FORMAT (44X,7HAZIMUTH/43X,9HDEVIATION,8X,9HELEVATION/ PRIN089
1 19X,16HHEIGHT RANGE,1X,2(5X,12HXMTR LOCAL),5X,26HPOLARIZATIPRIN090
2ON GROUP PATH,5A6,A5) PRIN091
write(2,1150) (UNIT(NN),NN=1,NP1) PRIN092
1150 FORMAT (13X,2(8X,2HKM),2X,2(6X,3HDEG,5X,3HDEG),6X,12HREAL IMAG,PRIN093
1 7X,2HKM,4X,3(4X,A6,2X)) PRIN094
IF (RAYSET.EQ.0.d0) GO TO 12 PRIN095
C********* PUNCH A TRANSMITTER RAYSET PRIN096
TLOND=TLON*DEGS PRIN097
IF (TLOND.LT.0.d0) TLOND=TLOND+360.d0 PRIN098
TLATD=TLAT*DEGS PRIN099
IF (TLATD.LT.0.d0) TLATD=TLATD+360.d0 PRIN100
AZ=AZ1*DEGS PRIN101
EL=BETA*DEGS PRIN102
NHOP=HOP PRIN103
ccccc write(7,1200)
c a ID(1),TYPE(NTYP),XMTRH,TLATD,TLOND,RCVRH,F,AZ,EL,POLAR PRIN104
c 1,NHOP,1HT PRIN105
1200 FORMAT (A3,A1,4PF9.0,2P2F6.0,4P2F9.0,5P2F10.0,5X,2P2F5.0,I1,A1) PRIN106
C********* PRIN107
12 V=0.d0 PRIN108
IF (N2.NE.0.d0) V=(R(4)**2+R(5)**2+R(6)**2)/N2-1.d0 PRIN109
H=R(1)-EARTHR PRIN110
STH=dSIN(THETA) PRIN111
CTH=dSIN(PID2-THETA) PRIN112
C********* CARTESIAN COORDINATES OF RAY POINT, ORIGIN AT TRANSMITTER PRIN113
XP=R(1)*STH*dSIN(PID2-PHI)-XR PRIN114
YP=R(1)*STH*dSIN(PHI)-YR PRIN115
ZP=R(1)*CTH-ZR PRIN116
C********* CARTESIAN COORDINATES OF RAY POINT, ORIGIN AT TRANSMITTER ANDPRIN117
C********* ROTATED PRIN118
EPS=XP*G1(1,1)+YP*G1(1,2)+ZP*G1(1,3) PRIN119
ETA=XP*G1(2,1)+YP*G1(2,2)+ZP*G1(2,3) PRIN120
ZETA=XP*G1(3,1)+YP*G1(3,2)+ZP*G1(3,3) PRIN121
RCE2=ETA**2+ZETA**2 PRIN122
RCE=dSQRT(RCE2) PRIN123
C********* GROUND RANGE PRIN124
RANGE=EARTHR*dATAN2(RCE,EARTHR+EPS+XMTRH) PRIN125
C********* ANGLE OF WAVE NORMAL WITH LOCAL HORIZONTAL PRIN126
ELL=dATAN2(R(4),dSQRT(R(5)**2+R(6)**2))*DEGS PRIN127
c
c see if ray was perpendicular to dipole field line between
c successive prints - check for sign change in cosine
ell_dip = ell/degs + dipole_dip
cos_ell = cos(ell_dip)
ccc print*,' dip, local elev ', dipole_dip*degs,ell
if(cos_ell_save / cos_ell .lt. 0.) print*,' Ray perp to field'
if(cos_ell_save / cos_ell .lt. 0.)
a write(2,2222) dipole_dip*degs,ell
2222 format(5x,'Ray was perp to the field',
a ' -- DIP & ELEV =',2f8.3)
cos_ell_save = cos_ell
c
C********* STRAIGHT LINE DISTANCE FROM TRANSMITTER TO RAY POINT PRIN128
SR=dSQRT(RCE2+EPS**2) PRIN129
IF (NP.LT.1) GO TO 16 PRIN130
DO 15 I=1,NP PRIN131
NN=NPR(I) PRIN132
15 RPRINT(I)=R(NN) PRIN133
16 IF (SR.GE.1.d-6) GO TO 20 PRIN134
C********* TOO CLOSE TO TRANSMITTER TO CALCULATE DIRECTION FROM PRIN135
C********* TRANSMITTER PRIN136
write(2,1500) V,NWHY,H,RANGE,ELL,POLAR,T,(RPRINT(NN),NN=1,NP1) PRIN137
1500 FORMAT (1X,1pE7.0,A8,0pF10.4,F11.4,26X,F8.3,F9.3,F8.3,8F12.4) PRIN138
GO TO 40 PRIN139
C********* ELEVATION ANGLE OF RAY POINT FROM TRANSMITTER PRIN140
20 EL=dATAN2(EPS,RCE)*DEGS PRIN141
IF (RCE.GE.1.d-6) GO TO 30 PRIN142
C********* NEARLY DIRECTLY ABOVE OR BELOW TRANSMITTER. CAN NOT CALCULATEPRIN143
C********* AZIMUTH DIRECTION FROM TRANSMITTER ACCURATELY PRIN144
write(2,2500) V,NWHY,H,RANGE,EL,ELL,POLAR,T,(RPRINT(NN),NN=1,NP1) PRIN145
2500 FORMAT (1X,1pE7.0,A8,0pF10.4,F11.4,17X,F9.3,F8.3,F9.3,F8.3, PRIN146
1 8F12.4) PRIN147
GO TO 40 PRIN148
C********* AZIMUTH ANGLE OF RAY POINT FROM TRANSMITTER PRIN149
30 ANGA=dATAN2(ETA,ZETA) PRIN150
AZDEV=180.d0-dMOD(540.d0-(AZ1-ANGA)*DEGS,360.d0) PRIN151
IF (R(5).NE.0.d0.OR.R(6).NE.0.d0) GO TO 34 PRIN152
C********* WAVE NORMAL IS VERTICAL, SO AZIMUTH DIRECTION CANNOT BE PRIN153
C********* CALCULATED PRIN154
write(2,3000)V,NWHY,H,RANGE,AZDEV,EL,ELL,POLAR,T,(RPRINT(NN),NN=1,PRIN155
1 NP1) PRIN156
3000 FORMAT (1X,1pE7.0,A8,0pF10.4,F11.4,F9.3,8X,F9.3,F8.3,F9.3,F8.3, PRIN157
1 8F12.4) PRIN158
GO TO 40 PRIN159
34 ANA=ANGA-ALPH PRIN160
SANA=dSIN(ANA) PRIN161
SPHI=SANA*STHR/STH PRIN162
CPHI=-dSIN(PID2-ANA)*dSIN(PID2-(PHI-PHIR))+SANA*dSIN(PHI-PHIR) PRIN163
1 *CTHR PRIN164
AZA=180.d0-dMOD (540.d0
a -(dATAN2(SPHI,CPHI)-dATAN2(R(6),R(5)))*DEGS,360.d0) PRIN165
cc type*,' PRINTR t ',t
write(2,3500)
a V,NWHY,H,RANGE,AZDEV,AZA,EL,ELL,POLAR,T,(RPRINT(NN),NN PRIN166
1 =1,NP1) PRIN167
3500 FORMAT (1X,1pE7.0,A8,0pF10.4,F11.4,2(F9.3,F8.3),F9.3,F8.3, PRIN168
1 8F12.4) PRIN169
C********* PRIN170
40 LINES=LINES+1 PRIN171
IF (NP.LE.3) GO TO 45 PRIN172
C********* ADDITIONAL LINE TO PRINT REMAINING DEPENDENT INTEGRATION PRIN173
C********* VARIABLES PRIN174
write(2,4000) (RPRINT(NN),NN=4,NP) PRIN175
4000 FORMAT (99X,3F12.4) PRIN176
LINES=LINES+1 PRIN177
45 IF (CARD.EQ.0.d0) RETURN PRIN178
C PRIN179
C********* PUNCH A RAYSET PRIN180
IF (AZDEV.LT.-90.d0) AZDEV=AZDEV+360.d0 PRIN181
IF (AZA.LT.-90.d0) AZA=AZA+360.d0 PRIN182
TDEV=T-SR PRIN183
NR=6 PRIN184
IF (W(57).EQ.0.d0) GO TO 47 PRIN185
C********* PHASE PATH PRIN186
NR=NR+1 PRIN187
PDEV=R(NR)-SR PRIN188
print*, ' Phase path ', pdev
47 IF (W(58).EQ.0.d0) GO TO 48 PRIN189
C********* ABSORPTION PRIN190
NR=NR+1 PRIN191
ABSORB=R(NR) PRIN192
C********* DOPPLER SHIFT PRIN193
48 IF (W(59).NE.0.d0) DOPP=R(NR+1) PRIN194
c write(7,4500)HPUNCH,RANGE,AZDEV,AZA,ELL,SR,TDEV,PDEV,ABSORB,DOPP, PRIN195
c 1 POLAR,IHOP,NWHY PRIN196
4500 FORMAT (4P2F9.0,3P3F6.0,3PF8.0,3P4F6.0,2P2F5.0,I1,A1) PRIN197
RETURN PRIN198
END PRIN199-