PROGRAM ESTPM (INPUT,OUTPUT,TAPE1=INPUT,TAPE2,TAPE3,TAPE4, . TAPE10=OUTPUT,TAPE11,TAPE12) COMMON A(1) COMPLEX CZ,ZZ,EE,FF,GG,CF,CZZ COMMON /REF/ X0,Y0,U0,V0,SX,SY,SU,SV,QP COMMON /GDL/ GW1,GW2,GN1,GN2,GRD,SGD,SCV,RCV,SBL DIMENSION NAME(8),FMT(8),BUFFER(1024) DATA (FMT(I),I=1,8) / "(A3,4X,A7,","1X,A10,A5,","10X,2(F2.0", . ",1X),F8.5,","1X,F3.0,1X",",F2.0,1X,F","8.5,F8.3,A", . "1) " / 900 FORMAT (8A10) AAA 7 901 FORMAT (I2,3I1,4F5.1,9F5.3,9X,A1) 902 FORMAT (1H1,/////20X,"EEEEEE SSSS TTTTTT PPPPP MM MM",25X, . 30("*"),/20X,"EE SS TT PP PP MMM MMM",25X, . "* J.A.R. BLAIS - JULY 1979 *",/20X,"EEEE SSS T", . "T PPPPP MM M MM",25X,"* PHYSICAL GEODESY SECTION *", . /20X,"EE SS TT PP MM MM",25X, . "* GEODETIC SURVEY OF CANADA *",/20X,"EEEEEE SSSS ", . "TT PP MM MM",25X,30("*")//) 905 FORMAT (//10X,"IDENTIFICATION OF PROJECT",10X,8A10/) 910 FORMAT (//10X,"ALTERNATE INPUT FORMAT",5X,8A10) 914 FORMAT (/10X,"THE GRID LATITUDE LIMITS ARE ",F6.1," AND ",F6.1, . " (+ NORTH)",//10X,"THE GRID LONGITUDE LIMITS ARE ",F6.1, . " AND ",F6.1," (+ WEST)") 915 FORMAT (/10X,"WRONG LIMITS FOR THE GRID") 921 FORMAT (///10X,"THE LOCAL ORIGIN IN LATITUDE IS ",F8.3, AAA 11 . " DEGREES (+ NORTH)") AAA 12 922 FORMAT (/10X,"THE LOCAL ORIGIN IN LONGITUDE IS ",F8.3, AAA 13 . " DEGREES (+ WEST)") 924 FORMAT (/10X,"THE NUMBER OF GEODETIC STATIONS USED IS ",I10) 925 FORMAT (//10X,"TRANSFORMATION FROM NAD27 TO ATS77 ",A10) 926 FORMAT (//10X,"TRANSFORMATION FROM ATS77 TO NAD27 ",A10) 931 FORMAT (//10X,"THE SCALE FACTOR IN LATITUDE IS ",F8.3) AAA 17 932 FORMAT (/10X,"THE SCALE FACTOR IN LONGITUDE IS ",F8.3) AAA 18 934 FORMAT (//10X,"THE NORMALIZED UNIT DISTANCE IN LATITUDE IS ", . F10.2,3X,"METRES") 935 FORMAT (//10X,"THE NORMALIZED UNIT DISTANCE IN LONGITUDE IS ", . F10.2,3X,"METRES") 940 FORMAT (///10X,"THE COMPUTED COEFFICIENTS ARE "/) AAA 20 942 FORMAT (20X,I5,5X,2F20.10) BBB 2 945 FORMAT (1H1///5X,"GEODETIC STATIONS",13X,"COORDINATES",13X, . "DIFFERENCES",22X,"ESTIMATES",15X,"RESIDUALS") 948 FORMAT (///10X,"THE DEGREE OF THE EXPANSION IS ",I5, AAA 22 . ///10X,"THE NUMBER OF COEFFICIENTS IS ",I6) AAA 23 949 FORMAT (///10X,"*** WARNING *** THE MAXIMUM NUMBER OF ", AAA 24 . "COEFFICIENTS IS ONLY ",I5) AAA 25 950 FORMAT (//5X,"NUMBER",5X,"NAME",15X,"LAT.",6X,"LONG.",8X, . "D(LAT)",6X,"D(LONG)",8X,"E(LAT)",10X,"E(LONG)",6X, . "R(LAT)",1X,"R(LONG)",2X,"LENGTH"/) 951 FORMAT (34X,"(DEG)",6X,"(DEG)",8X,"(SEC)",7X,"(SEC)",4X, . 3(4X,"(SEC)",2X,"(SEC)"),6X,"(M)"/) 952 FORMAT (4X,2(1X,A10),A5,1X,2(F8.3,2X),3X,2(F8.3,4X),3(2X,2F7.3), . 3X,F7.3) 955 FORMAT (//10X,"THE RMS VALUE IN LATITUDE IS ",F10.3," SECONDS") 956 FORMAT (//10X,"THE RMS VALUE IN LONGITUDE IS ",F10.3," SECONDS") 960 FORMAT (1H1,///5X,"IDENTIFICATION OF PROJECT",10X,8A10,/////5X, . "TRANSFORMATION OF OTHER POINTS",46X,"APPLIED CORRECTIONS", . " IN SECONDS") 961 FORMAT (//5X,"NUMBER",5X,"NAME",15X,"LATITUDE",13X,"LONGITUDE", . 19X,"LATITUDE",8X,"LONGITUDE"/) 962 FORMAT (5X,A7,1X,A10,A5,2I5,F10.5,3X,2I5,F10.5,10X,2(2F7.3,2X)) 965 FORMAT (A3,4X,A7,1X,A10,A5,9X,2I3,F9.5,1X,2I3,1X,F8.5,F8.3,A1) 970 FORMAT (///10X,"THE REQUESTED CENTRAL MEMORY IS ",O6," OCTAL ", AAA 43 . "WORDS") AAA 44 990 FORMAT (/////10X,"THE REQUESTED CENTRAL MEMORY IS TOO SMALL BY ", AAA 45 . I10," DECIMAL WORDS") AAA 46 X0=Y0=U0=V0=SX=SY=SU=SV=0.0 ICR = 1 AAA 48 IT2 = 2 AAA 49 IT3 = 3 IT4 = 4 IPR = 10 ITP = 11 ITQ = 12 SGN = - 1.0 AAA 52 CODE = 3H 4 AAA 53 QODE = 3H 04 SCRD = ASIN(1.0)/90.0 CALL JCMEM (IPG,ICOM,JCM) CALL KOREA (JCM,ICOM) CALL COMSAVE (JCM) C*****IDENTIFICATION OF THE PROJECT AAA 54 READ (ICR,900) (NAME(I),I=1,8) AAA 55 WRITE (IPR,902) WRITE (IPR,905) (NAME(I),I=1,8) READ (ICR,901) NS,K1,K2,K3,GW1,GW2,GN1,GN2,SGD,SCV,GRD,RCV, . SBN,SBS,SBI,RUI,RUE,KF IF (KF.EQ.1H ) GO TO 1 READ (ICR,900) (FMT(I),I=1,8) WRITE (IPR,910) (FMT(I),I=1,8) 1 IF (NS.EQ.99) READ (ITP) NR,NP,K1,RUI,X0,Y0,U0,V0,SX,SY,SU,SV IF (SGD.LE.0.0) SGD = 1.0 IF (SCV.LE.0.0) SCV = 1.0 IF (GRD.LE.0.0) GRD = 1.0 IF (SBN.LE.0.0) SBN = 0.21 IF (SBS.LE.0.0) SBS = 0.14 IF (SBI.LE.0.0) SBI = 0.07 IF (RCV.LE.0.0) RCV = 3.0 IF (RUI.EQ.0.0) RUI = 0.1 IF (RUI.LT.0.0) RUI = 0.0 IF (RUE.LE.0.0) RUE = 0.1 SBL = 2.0*SBN/3.0 IF (GW1.NE.GW2) GO TO 2 IF (GN1.NE.GN2) GO TO 2 WRITE(IPR,914) GN1,GN2,GW1,GW2 WRITE (IPR,915) STOP 2 IF (GW1.GT.GW2) GO TO 3 GW0 = GW1 GW1 = GW2 GW2 = GW0 3 IF (GN1.LT.GN2) GO TO 4 GN0 = GN1 GN1 = GN2 GN2 = GN0 4 IF (NS.NE.99) NR = NS NF = NR + 1 NC = 2*NF BBB 5 WRITE (IPR,948) NR,NC BBB 6 WRITE (IPR,970) JCM KCM = JCM - LOCF(A(1)) CKT = 10HNO IF (K1.EQ.1) CKT = 10HYES WRITE (IPR,925) CKT CKT = 10HNO IF (K1.EQ.2) CKT = 10HYES WRITE (IPR,926) CKT IF (NS.EQ.99) GO TO 10 C*****LATITUDE IN DEGREES (+ NORTH) AAA 74 C*****LONGITUDE IN DEGREES (+ WEST) ITS = IT2 5 READ (ITS,FMT) D0,D1,D2,D3,U1,U2,U3,V1,V2,V3,HH,HC IF (EOF(ITS).NE.0.0) GO TO 8 IF (KF.EQ.1HZ) GO TO 6 IF (D0.EQ.3H 4.OR.D0.EQ.3H 04) GO TO 6 IF (D0.EQ.3H 5.OR.D0.EQ.3H 05) GO TO 6 IF (D0.EQ.3H 6.OR.D0.EQ.3H 06) GO TO 6 IF (D0.EQ.3H 96) GO TO 6 GO TO 5 6 IF (ITS.EQ.2) GO TO 7 IF (K1.EQ.2) CALL NADATS (U1,U2,U3,V1,V2,V3,HH) IF (K1.EQ.3) CALL ATSNAD (U1,U2,U3,V1,V2,V3,HH) 7 X = SGN*(U1 + U2/60.0 + U3/3600.0) Y = SGN*(V1 + V2/60.0 + V3/3600.0) IF (ABS(X).LT.GN1) GO TO 5 IF (ABS(X).GT.GN2) GO TO 5 IF (ABS(Y).GT.GW1) GO TO 5 IF (ABS(Y).LT.GW2) GO TO 5 WRITE (ITP) D1,D2,D3,X,Y GO TO 5 8 IF (ITS.EQ.IT3) GO TO 9 ITS = IT3 SGN = 1.0 GO TO 5 9 REWIND ITP CALL DATAST (ITP,ITQ,JCM) NP = INT (QP+0.1) C*****LOCAL ORIGIN AS CENTRE OF GRAVITY OF DATA AAA 88 10 WRITE (IPR,921) X0 WRITE (IPR,922) Y0 AAA 93 WRITE (IPR,914) GN1,GN2,GW1,GW2 WRITE (IPR,924) NP IF (NP.EQ.0) STOP LCM = NF*(NF+1) + 4*NF + 5*NP - KCM IF (LCM.GT.0) WRITE (IPR,990) LCM BBB 11 IF (LCM.GT.0) STOP BBB 12 N1 = 1 BBB 13 N2 = 2*NP + 1 BBB 14 N3 = 4*NP + 1 BBB 15 N4 = 5*NP + 1 BBB 16 N5 = N4 + 2*NF N6 = N5 + 2*NF NM = (NF*(NF+1))/2 IF (NS.EQ.99) GO TO 11 CALL LSQSLN (A(N1),A(N2),A(N3),A(N4),A(N5),A(N6),NP,NR,NF,NM,ITP) C*****SCALE FACTORS FROM THE DISPERSION OF THE DATA AAA 102 11 WRITE (IPR,931) SX WRITE (IPR,932) SY AAA 104 RNX = RUI*111319.456/SX RNY = RUI*111319.456*COS(X0*SCRD)/SY WRITE (IPR,934) RNX WRITE (IPR,935) RNY IF (K2.EQ.0) GO TO 12 DR = FLOAT (NR) SGN = - SGD/COS(GN1*SCRD) ASW = (GW2-GW1)/SGN ASN = (GN2-GN1)/SGD CALL PLOTS (BUFFER,1024) CALL PLOT (0.0,-30.0,-3) CALL PLOT (4.0,2.0,-3) CALL SYMBOL (0.0,-2.0,SBN,NAME,0.0,80) CALL SYMBOL (0.0,-1.0,SBL,"DEGREE OF ESTIMATION = ",0.0,23) CALL NUMBER (5.0,-1.0,SBL,DR,0.0,-1) CALL SYMBOL (0.0,-1.5,SBL,"VECTOR SCALE IN SECONDS PER INCH = ", . 0.0,35) CALL NUMBER (5.0,-1.5,SBL,SCV,0.0,1) CALL AXIS (0.0,0.0,27HLONGITUDE POSITIVE WESTWARD,-27,ASW,0.0, . GW1,SGN) CALL AXIS (0.0,0.0,27HLATITUDE POSITIVE NORTHWARD,27,ASN,90.0, . GN1,SGD) CALL AXIS (ASW,0.0,27HLATITUDE POSITIVE NORTHWARD,-27,ASN,90.0, . GN1,SGD) SGN = - SGD / COS(GN2*SCRD) GWN = (GW1+GW2-ASW*SGN) / 2.0 CALL AXIS (0.0,ASN,27HLONGITUDE POSITIVE WESTWARD,27,ASW,0.0, . GWN,SGN) 12 IF (NS.NE.99) REWIND ITP WRITE (IPR,940) AAA 108 IF (NS.EQ.99) READ (ITP) (A(I),I=1,N5) DO 14 I=1,NF II = I - 1 BBB 22 C1 = A(N4+2*I-2) C2 = A(N4+2*I-1) 14 WRITE (IPR,942) II,C1,C2 RP = RQ = PU = PV = QU = QV = 0.0 BBB 24 IF (NS.EQ.99) GO TO 21 WRITE (IPR,945) WRITE (IPR,950) AAA 121 WRITE (IPR,951) AAA 122 DO 20 I=1,NP READ (ITP) PT,PT1,PT2,X,Y,U,V,W CZ = CMPLX (SX*(X-X0),SY*(Y-Y0)) FF = CMPLX (SU*(U-U0)-A(N2+2*I-2),SV*(V-V0)-A(N2+2*I-1)) QQ = 0.0 GG = CMPLX (0.0,0.0) IF (RUI.EQ.0.0) GO TO 18 DO 15 J=1,NP ZZ = CMPLX (A(N1+2*J-2),A(N1+2*J-1)) EE = CMPLX (A(N2+2*J-2),A(N2+2*J-1)) DD = (CZ-ZZ)*CONJG(CZ-ZZ) / (RUI*RUI) IF (DD.GT.675.0) DD = 675.0 WW = A(N3+J-1) QW = WW*EXP(-DD) GG = (QQ*GG + QW*EE) / (QQ + QW) 15 QQ = QQ + QW 18 FU = REAL (FF) / SU + U0 FV = AIMAG (FF) / SV + V0 GU = REAL (GG) / SU GV = AIMAG (GG) / SV WU = WV = 1.0 RU = RV = 0.0 BBB 27 RU = FU + GU - U QU = (PU*QU + WU*RU*RU) / (PU + WU) AAA 134 PU = PU + WU AAA 135 RP = RP + RU BBB 29 RV = FV + GV - V QV = (PV*QV + WV*RV*RV) / (PV + WV) AAA 138 PV = PV + WV AAA 139 RQ = RQ + RV BBB 30 RR = SQRT (RU*RU + (RV*COS(X*SCRD))**2)*30.92209533 WRITE (IPR,952) PT,PT1,PT2,X,Y,U,V,FU,GU,FV,GV,RU,RV,RR IF (K2.EQ.0) GO TO 20 SGN = - SGD / COS(X*SCRD) GWN = (GW1+GW2-ASW*SGN) / 2.0 PW = (Y-GWN) / SGN PN = (X-GN1) / SGD QW = PW - V*COS(X*SCRD)/SCV QN = PN + U / SCV IF (K2.EQ.2) GO TO 19 CALL SYMBOL (PW+SBI,PN-SBI,SBI,PT,0.0,10) PQ = (PW-QW)**2 + (PN-QN)**2 IF (PQ.LT.RCV*RCV) GO TO 19 CALL SYMBOL (PW,PN,SBS,14,0.0,-1) GO TO 20 19 CALL SYMBOL (PW,PN,SBS,2,0.0,-1) CALL PLOT (PW,PN,3) CALL PLOT (QW,QN,2) 20 CONTINUE 21 CONTINUE IF (K3.EQ.1) CALL GRID (A(N1),A(N2),A(N3),A(N4),NP,NF,RUI,RUE) IF (K3.EQ.2) CALL GRID (A(N1),A(N2),A(N3),A(N4),NP,NF,RUI,-1.1) IF (NS.EQ.99) GO TO 48 REWIND IT2 QU = SQRT (QU) QV = SQRT (QV) WRITE (IPR,955) QU AAA 146 WRITE (IPR,956) QV AAA 147 48 KPT = 0 QQQ = 0.0000001 50 READ (IT4,FMT) D0,D1,D2,D3,U1,U2,U3,V1,V2,V3,HH,HC IF (EOF(IT4).NE.0.0) GO TO 90 IF (KF.EQ.1HZ) GO TO 51 IF (D0.EQ.3H 4.OR.D0.EQ.3H 04) GO TO 51 IF (D0.EQ.3H 5.OR.D0.EQ.3H 05) GO TO 51 IF (D0.EQ.3H 6.OR.D0.EQ.3H 06) GO TO 51 IF (D0.EQ.3H 96) GO TO 51 GO TO 50 51 IF (K1.EQ.2) CALL NADATS (U1,U2,U3,V1,V2,V3,HH) IF (K1.EQ.3) CALL ATSNAD (U1,U2,U3,V1,V2,V3,HH) X = U1 + U2/60.0 + U3/3600.0 Y = V1 + V2/60.0 + V3/3600.0 IF (ABS(X).LT.GN1) GO TO 50 IF (ABS(X).GT.GN2) GO TO 50 IF (ABS(Y).GT.GW1) GO TO 50 IF (ABS(Y).LT.GW2) GO TO 50 CZ = CMPLX (SX*(X-X0),SY*(Y-Y0)) IF (CABS(CZ).GT.(1.0+RUE)) GO TO 50 KPT = KPT + 1 IF (KPT.GT.1) GO TO 52 WRITE (IPR,960) (NAME(I),I=1,8) WRITE (IPR,961) 52 CF = CMPLX (A(N4),A(N4+1)) CZZ = CMPLX (1.0,0.0) DO 54 K=2,NF CZZ = CZ*CZZ 54 CF = CF + CZZ*CMPLX (A(N4+2*K-2),A(N4+2*K-1)) QQ = 0.0 CZZ = CMPLX (0.0,0.0) IF (RUI.EQ.0.0) GO TO 61 DO 60 I=1,NP ZZ = CMPLX (A(N1+2*I-2),A(N1+2*I-1)) EE = CMPLX (A(N2+2*I-2),A(N2+2*I-1)) DD = (CZ-ZZ)*CONJG(CZ-ZZ) / (RUI*RUI) IF (DD.GT.675.0) DD = 675.0 WW = A(N3+I-1) QW = WW*EXP(-DD) CZZ = (QQ*CZZ + QW*EE) / (QQ + QW) 60 QQ = QQ + QW 61 FU = REAL(CF) / SU + U0 FV = AIMAG (CF) / SV + V0 GU = REAL (CZZ) / SU GV = AIMAG (CZZ) / SV XX = 3600.0*X + FU + GU YY = 3600.0*Y + FV + GV IU1 = IFIX ((XX+QQQ)/3600.0) IV1 = IFIX ((YY+QQQ)/3600.0) IU2 = IFIX ((XX+QQQ)/60.0) - IU1*60 IV2 = IFIX ((YY+QQQ)/60.0) - IV1*60 U3 = XX - 60*IU2 - 3600*IU1 V3 = YY - 60*IV2 - 3600*IV1 WRITE (IPR,962) D1,D2,D3,IU1,IU2,U3,IV1,IV2,V3,FU,GU,FV,GV WRITE (IT2,965) D0,D1,D2,D3,IU1,IU2,U3,IV1,IV2,V3,HH,HC IF (K2.EQ.0) GO TO 50 SGN = - SGD / COS(X*SCRD) GWN = (GW1+GW2-ASW*SGN) / 2.0 PW = (Y-GWN) / SGN PN = (X-GN1) / SGD QW = PW - (FV+GV)*COS(X*SCRD)/SCV QN = PN + (FU+GU)/SCV IF (K2.EQ.2) GO TO 62 CALL SYMBOL (PW+SBI,PN-SBI,SBI,D1,0.0,10) PQ = (PW-QW)**2 + (PN-QN)**2 IF (PQ.LT.RCV*RCV) GO TO 62 CALL SYMBOL (PW,PN,SBS,10,0.0,-1) GO TO 50 62 CALL SYMBOL (PW,PN,SBS,1,0.0,-1) CALL PLOT (PW,PN,3) CALL PLOT (QW,QN,2) GO TO 50 90 IF (K2.NE.0) CALL PLOT (15.0,0.0,999) IF (NS.EQ.99) GO TO 100 REWIND ITP WRITE (ITP) NR,NP,K1,RUI,X0,Y0,U0,V0,SX,SY,SU,SV WRITE (ITP) (A(I),I=1,N5) 100 REWIND IT2 REWIND IT3 REWIND ITP REWIND ITQ CALL COMREST STOP END AAA 149 SUBROUTINE NADATS (U1,U2,U3,V1,V2,V3,HE) C*****TRANSFORMATION FROM NAD27 TO ATS77 C*****INPUT AND OUTPUT LATITUDES ARE POSITIVE NORTHWARD C*****INPUT AND OUTPUT LONGITUDES ARE POSITIVE WESTWARD C*****U1,U2,U3 = LATITUDE IN DEG., MIN., SEC. C*****V1,V2,V3 = LONGITUDE IN DEG., MIN., SEC. C*****HE = ELLIPSOIDAL HEIGHT IN METRES C*****REF. HEISKANEN AND MORITZ = PHYSICAL GEODESY P.207, EQN.5-55 SCRD = ASIN(1.0) / 90.0 AAA = 6378135.0 FFF = 1.0 / 298.257 AA = 6378206.4 BB = 6356583.8 FF = (AA-BB)/AA DA = AAA - AA DF = FFF - FF DX = 15.0 DY = - 165.0 DZ = - 175.0 UU = (U1 + U2/60.0 + U3/3600.0)*SCRD VV = 4.0*ASIN(1.0) - (V1 + V2/60.0 + V3/3600.0)*SCRD SU = SIN (UU) SV = SIN (VV) CU = COS (UU) CV = COS (VV) DU = (SU*CV*DX + SU*SV*DY - CU*DZ + 2.0*AA*SU*CU*DF) / AA DV = (SV*DX - CV*DY) / (AA*CU) DH = - CU*CV*DX - CU*SV*DY - SU*DZ - DA + AA*SU*SU*DF UU = UU + DU VV = 4.0*ASIN(1.0) - VV - DV U1 = AINT(UU/SCRD) U2 = AINT(UU*60.0/SCRD) - 60.0*U1 U3 = UU*3600.0/SCRD - 3600.0*U1 - 60.0*U2 V1 = AINT(VV/SCRD) V2 = AINT(VV*60.0/SCRD) - 60.0*V1 V3 = VV*3600.0/SCRD - 3600.0*V1 - 60.0*V2 HE = HE + DH RETURN END SUBROUTINE ATSNAD (U1,U2,U3,V1,V2,V3,HE) C*****TRANSFORMATION FROM ATS77 TO NAD27 C*****INPUT AND OUTPUT LATITUDES ARE POSITIVE NORTHWARD C*****INPUT AND OUTPUT LONGITUDES ARE POSITIVE WESTWARD C*****U1,U2,U3 = LATITUDE IN DEG., MIN., SEC. C*****V1,V2,V3 = LONGITUDE IN DEG., MIN., SEC. C*****HE = ELLIPSOIDAL HEIGHT IN METRES C*****REF. HEISKANEN AND MORITZ = PHYSICAL GEODESY P.207 SCRD = ASIN(1.0) / 90.0 AAA = 6378135.0 FFF = 1.0 / 298.257 AA = 6378206.4 BB = 6356583.8 FF = (AA-BB)/AA DA = AA - AAA DF = FF - FFF DX = - 15.0 DY = 165.0 DZ = 175.0 UU = (U1 + U2/60.0 + U3/3600.0)*SCRD VV = 4.0*ASIN(1.0) - (V1 + V2/60.0 + V3/3600.0)*SCRD SU = SIN (UU) SV = SIN (VV) CU = COS (UU) CV = COS (VV) DU = (SU*CV*DX + SU*SV*DY - CU*DZ + 2.0*AAA*SU*CU*DF) / AAA DV = (SV*DX - CV*DY) / (AAA*CU) DH = - CU*CV*DX - CU*SV*DY - SU*DZ - DA + AAA*SU*SU*DF UU = UU + DU VV = 4.0*ASIN(1.0) - VV - DV U1 = AINT(UU/SCRD) U2 = AINT(UU*60.0/SCRD) - 60.0*U1 U3 = UU*3600.0/SCRD - 3600.0*U1 - 60.0*U2 V1 = AINT(VV/SCRD) V2 = AINT(VV*60.0/SCRD) - 60.0*V1 V3 = VV*3600.0/SCRD - 3600.0*V1 - 60.0*V2 HE = HE + DH RETURN END SUBROUTINE DATAST (ITP,ITQ,JCM) COMMON /REF/ X0,Y0,U0,V0,SX,SY,SU,SV,QP COMMON /SPAN/ XA,XB,YA,YB CALL COMREST CALL SMSORT (50) CALL SMFILE ("SORT","BINARY",ITP,"REWIND") CALL SMFILE ("OUTPUT","BINARY",ITQ,"REWIND") CALL SMKEY (1,1,10,0,"LOGICAL") CALL SMEND CALL COMSAVE (JCM) QP = 0.0 W = 1.0 QT = BLK = 10H 1 READ (ITQ) PT,PT1,PT2,X,Y IF (EOF(ITQ).NE.0.0) PT = BLK IF (PT.EQ.BLK) GO TO 10 IF (PT.NE.QT) GO TO 10 2 IF (X.GT.0.0) GO TO 5 X1 = (Q1*X1+X) / (Q1+W) Y1 = (Q1*Y1+Y) / (Q1+W) Q1 = Q1 + W GO TO 1 5 X2 = (Q2*X2+X) / (Q2+W) Y2 = (Q2*Y2+Y) / (Q2+W) Q2 = Q2 + W GO TO 1 10 IF (QT.EQ.BLK) GO TO 12 IF (Q1*Q2.LT.0.5) GO TO 12 X1 = - X1 Y1 = - Y1 U = 3600.0*(X2-X1) V = 3600.0*(Y2-Y1) WRITE (ITP) QT,QT1,QT2,X1,Y1,U,V,W QQ = QP + W X0 = (QP*X0+X1) / QQ Y0 = (QP*Y0+Y1) / QQ U0 = (QP*U0+U) / QQ V0 = (QP*V0+V) / QQ QP = QQ IF (QP.GT.1.5) GO TO 20 XA = XB = X1 YA = YB = Y1 20 XA = AMIN1 (XA,X1) XB = AMAX1 (XB,X1) YA = AMIN1 (YA,Y1) YB = AMAX1 (YB,Y1) SU = AMAX1 (SU,U*U) SV = AMAX1 (SV,V*V) 12 IF (PT.EQ.BLK) GO TO 40 QT = PT QT1 = PT1 QT2 = PT2 X1 = Y1 = Q1 = 0.0 X2 = Y2 = Q2 = 0.0 GO TO 2 40 IF (QP.GT.0.5) GO TO 50 SX=SY=SU=SV=0.0 ENDFILE ITP RETURN 50 IF (XA.NE.XB) GO TO 52 SX = 1.0 GO TO 53 52 SX = 1.0 / AMAX1(XB-X0,X0-XA) 53 IF (YA.NE.YB) GO TO 54 SY = 1.0 GO TO 55 54 SY = 1.0 / AMAX1(YB-Y0,Y0-YA) 55 IF (SU.GT.U0*U0) GO TO 56 SU = 1.0 GO TO 57 56 SU = 1.0 / SQRT(SU-U0*U0) 57 IF (SV.GT.V0*V0) GO TO 58 SV = 1.0 GO TO 59 58 SV = 1.0 / SQRT(SV-V0*V0) 59 REWIND ITP REWIND ITP RETURN AAA 222 END AAA 223 SUBROUTINE LSQSLN (ZZ,EE,WW,CC,DD,AA,NP,NR,NF,NM,ITP) COMMON /REF/ X0,Y0,U0,V0,SX,SY,SU,SV,QP COMMON /SPAN/ X1,X2,Y1,Y2 COMPLEX ZZ(NP),EE(NP),CC(NF),DD(NF),AA(NM),CF DIMENSION WW(NP) CCC 22 DO 1 I=1,NF CCC 23 1 DD(I) = CMPLX (0.0,0.0) CCC 24 DO 2 J=1,NM CCC 25 2 AA(J) = CMPLX (0.0,0.0) CCC 26 DO 5 I=1,NP CCC 27 READ (ITP) PT,PT1,PT2,X,Y,U,V,W ZZ(I) = CMPLX (SX*(X-X0),SY*(Y-Y0)) EE(I) = CMPLX (SU*(U-U0),SV*(V-V0)) WW(I) = W CC(1) = CMPLX (1.0,0.0) DD(1) = DD(1) + WW(I)*EE(I) CCC 29 AA(1) = AA(1) + WW(I) CCC 30 DO 3 J=2,NF CCC 31 J1 = KFN(1,J,NF) CCC 32 CC(J) = ZZ(I)*CC(J-1) 3 AA(J1) = AA(J1) + WW(I)*CC(J) IF (NF.EQ.1) GO TO 5 CCC 34 DO 4 K=2,NF CCC 35 DD(K) = DD(K) + WW(I)*EE(I)*CONJG(CC(K)) DO 4 L=K,NF CCC 37 KL = KFN(K,L,NF) CCC 38 4 AA(KL) = AA(KL) + WW(I)*CC(L)*CONJG(CC(K)) 5 CONTINUE CCC 40 CALL CXSLN (AA,DD,CC,NM,NF,NRT) CCC 41 DO 10 I=1,NP CF = CC(1) IF (NR.EQ.0) GO TO 10 DD(1) = CMPLX (1.0,0.0) DO 8 J=2,NF DD(J) = ZZ(I)*DD(J-1) 8 CF = CF + CC(J)*DD(J) 10 EE(I) = EE(I) - CF RETURN END CCC 43 SUBROUTINE CXSLN (AA,AB,CT,NM,NT,NRT) CCC 44 COMPLEX AA(NM),AB(NT),CT(NT),PP,QQ CCC 45 C*****FORWARD REDUCTION CCC 46 QQ = CSQRT (AA(1)) CCC 47 DO 1 I=1,NT CCC 48 1 AA(I) = AA(I) / QQ CCC 49 AB(1) = AB(1) / QQ CCC 50 I1 = 0 CCC 51 IF (NT.EQ.1) GO TO 7 CCC 52 DO 6 I=2,NT CCC 53 I1 = I - 1 CCC 54 I2 = I + 1 CCC 55 PP = CMPLX (0.0,0.0) CCC 56 DO 2 L=1,NT CCC 57 2 CT(L) = CMPLX (0.0,0.0) CCC 58 DO 4 K=1,I1 CCC 59 KI = KFN(K,I,NT) CCC 60 DO 3 J=I,NT CCC 61 KJ = KFN(K,J,NT) CCC 62 3 CT(J) = CT(J) + CONJG(AA(KI))*AA(KJ) CCC 63 4 PP = PP + CONJG(AA(KI))*AB(K) CCC 64 II = KFN(I,I,NT) CCC 65 IF (REAL(AA(II)-CT(I)).LE.0.001) GO TO 7 CCC 66 QQ = CSQRT (AA(II)-CT(I)) CCC 67 AA(II) = QQ CCC 68 IF (I.EQ.NT) GO TO 6 CCC 69 DO 5 M=I2,NT CCC 70 IM = KFN(I,M,NT) CCC 71 5 AA(IM) = AA(IM)/QQ - CT(M)/QQ CCC 72 6 AB(I) = AB(I)/QQ - PP/QQ CCC 73 C*****BACKWARD SUBSTITUTION CCC 74 7 DO 8 I=1,NT CCC 75 8 CT(I) = CMPLX (0.0,0.0) CCC 76 IF (I1.EQ.(NT-1)) I1 = NT CCC 77 NRT = I1 CCC 78 II = KFN(I1,I1,NT) CCC 79 CT(I1) = AB(I1) / AA(II) CCC 80 IF (I1.EQ.1) RETURN CCC 81 I2 = I1 - 1 CCC 82 DO 10 I=1,I2 CCC 83 J = I2 - I + 1 CCC 84 JJ = KFN(J,J,NT) CCC 85 CT(J) = AB(J) / AA(JJ) CCC 86 JP1 = J + 1 CCC 87 DO 9 K=JP1,I1 CCC 88 JK = KFN(J,K,NT) CCC 89 9 CT(J) = CT(J) - AA(JK)*CT(K)/AA(JJ) CCC 90 10 CONTINUE CCC 91 RETURN CCC 92 END CCC 93 SUBROUTINE GRID (ZZ,EE,WW,CFN,NP,NF,RUI,RUE) COMMON /REF/ X0,Y0,U0,V0,SX,SY,SU,SV COMMON /GDL/ GW1,GW2,GN1,GN2,GRD,SGD,SCV,RCV,SBL COMPLEX ZZ(NP),EE(NP),CFN(NF),CZ,CZZ,CF DIMENSION WW(NP) SCRD = ASIN(1.0) / 90.0 IF (GW1-GRD.LT.GW2) RETURN GN = GN1 + GRD IF (GN.GT.GN2) RETURN SGN = - SGD / COS(GN1*SCRD) ASW = (GW2-GW1) / SGN 1 GW = GW1 + GRD 2 GW = GW - GRD SGN = - SGD / COS(GN*SCRD) GWN = (GW1+GW2-ASW*SGN) / 2.0 PW = (GW - GWN) / SGN PN = (GN - GN1) / SGD CZ = CMPLX (SX*(GN-X0),SY*(GW-Y0)) CF = CFN(1) CZZ = CMPLX (1.0,0.0) DO 3 I=2,NF CZZ = CZ*CZZ 3 CF = CF + CFN(I)*CZZ QQ = 0.0 CZZ = CMPLX (0.0,0.0) IF (RUI.EQ.0.0) GO TO 6 DO 5 I=1,NP DD = (CZ-ZZ(I))*CONJG(CZ-ZZ(I)) / (RUI*RUI) IF (DD.GT.675.0) DD = 675.0 QW = WW(I)*EXP(-DD) CZZ = (QQ*CZZ + QW*EE(I)) / (QQ + QW) 5 QQ = QQ + QW 6 CF = CF + CZZ EU = REAL(CF)/SU + U0 EV = AIMAG(CF)/SV + V0 IF (CABS(CZ).GT.(1.0+RUE)) EU = EV = 0.0 QW = PW - EV*COS(GN*SCRD)/SCV QN = PN + EU/SCV PQ = (PW-QW)**2 + (PN-QN)**2 IF (PQ.LT.RCV*RCV) GO TO 8 CALL SYMBOL(PW,PN,SBL,11,0.0,-1) GO TO 10 8 CALL SYMBOL(PW,PN,SBL,3,0.0,-1) CALL PLOT (PW,PN,3) CALL PLOT (QW,QN,2) 10 IF (GW-GRD.GE.GW2) GO TO 2 GN = GN + GRD IF (GN.LT.GN2) GO TO 1 RETURN END FUNCTION KFN (I,J,N) AAA 462 KFN = (I-1)*N - (I*I-I)/2 + J AAA 463 RETURN AAA 464 END AAA 465 IDENT JCMEM ENTRY JCMEM USE // A BSS 1 SET A AS FIRST WORD IN BLANK COMMON USE 0 STATUS VFD 30/-1,29/0 SET -1 IN BITS 30-59 REST 0 JCMEM DATA 0 SA3 A1 SAVE ADDRESS OF F.P. 1 MEMORY CM,STATUS,R SA2 STATUS GET THE RETURNED CM VALUE IN BITS 30-59 AX2 30 SHIFT TO STANDARD INTEGER FORMAT SX7 A PUT ADDRESS OF BLANK COMMON IN F.P. 1 SA7 X3 STORE VALUE SA3 A3+1 SET ADDRESS TO F.P. 2 IX7 X2-X7 CALCULATE SIZE OF BLANK COMMON SX7 X7-1 SA7 X3 STORE IN F.P. 2 SA3 A3+1 SET ADDRESS OF F.P. 3 BX7 X2 STORE MAX FL IN F.P. 3 SA7 X3 EQ JCMEM END IDENT KOREA 000010 ENTRY KOREA 000020 STATUS BSS 1 000030 USE // 000040 A BSS 0 000050 USE 0 000060 KOREA DATA 0 000070 SA3 X1 000080 BX6 X3 000090 LX6 30 000100 SA6 STATUS 000110 SA3 A1 000120 MEMORY CM,STATUS,R,,1 SA2 STATUS 000140 AX2 30 000150 BX7 X2 000160 SA7 X3 STORE FIELD LENGTH 000170 SX5 A 000180 IX7 X7-X5 000190 SX7 X7-1 000200 SA3 A3+1 000210 SA7 X3 STORE KOMMON SIZE 000220 EQ KOREA 000230 END 000240 IDENT COMSAVE ENTRY COMSAVE ENTRY COMREST COMSAVE BSSZ 1 SAVE DABA TO FL SB5 B0 0 MEANS NO FL ARG ZR X1,NOFL IF NO ARG SA2 X1 FETCH FL VALUE SB5 X2 SAVE DESIRED FL NOFL BSS 0 SA1 65B RA+65B SX6 X1 SA6 CMMACTV PL X6,COMSAVE IF CMM NOT ACTIVE BX6 -X6 SA6 DABA SA6 OUT SA6 FIRST SET FET POINTERS FOR WRITE SA2 X6 GET DYNAMIC AREA HEADER SX7 X2 FIELD LENGTH AS KNOWN BY CMM SA7 FL SA7 LIMIT SX7 X7-1 FL-1 SA7 IN SA3 X7 GET WORD AT FL-1 BX7 X3 SAVE IT (CANT BE WRITTEN) SA7 FLSAVE SX2 24B RJ CIO WRITE DABA TO FL TO SAVE FILE SX2 50B RJ CIO REWIND SAVE FILE EQ B5,B0,COMSAVE RETURN IF NO FL ARG SX2 B5 DESIRED FL RJ CMEM ISSUE CMM MEM REQUEST FOR FL EQ COMSAVE RETURN SPACE 3 COMREST BSSZ 1 RESTORE DABA TO FL SA2 CMMACTV PL X2,COMREST RETURN IF NOTHING SAVED MX6 0 SA6 A2 INDICATE SAVE FILE IS OUT-OF-DATE SA2 FL CMM FL RJ CMEM ADJUST FL TO WHAT CMM WANTS SA2 FIRST SX6 X2 SA6 IN RESET FET POINTERS SA6 OUT SX2 10B RJ CIO READ SAVE FILE TO DABA TO FL SX2 50B RJ CIO REWIND SAVE FILE SA2 FLSAVE BX6 X2 SA3 FL SA6 X3-1 RESTORE WORD AT FL-1 EQ COMREST RETURN SPACE 3 CIO BSSZ 1 ISSUE CIO FUNCTION IN X2 SA3 FET MX0 42 BX6 X0*X3 CLEAR OLD FUNCTION BX6 X6+X2 OR IN NEW FUNCTION SA6 A3 TO FET SYSTEM CIO,RECALL,FET CALL CIO FOR FUNCTION EQ CIO RETURN SPACE 3 CMEM BSSZ 1 ISSUE CMM MEMORY REQUEST BX6 X2 LX6 30 SX2 4 BX6 X6+X2 FL+CMM FLAG SA6 MEMREQ MEMORY CM,MEMREQ,RECALL EQ CMEM RETURN SPACE 3 FET VFD 42/7LZZZZZUF,18/0 FIRST BSSZ 1 IN BSSZ 1 OUT BSSZ 1 LIMIT BSSZ 1 SPACE 3 FL BSSZ 1 FIELD LENGTH KNOWN TO CMM DABA BSSZ 1 DYNAMIC AREA BASE ADDRESS FLSAVE BSSZ 1 WORD AT FL-1 SAVED MEMREQ BSSZ 1 ARG FOR MEM CMMACTV BSSZ 1 NEGATIVE IF CMM ACTIVE (DABA TO FL SA END