C GALS MAINLINE TO ADJUST DIMENSIONS 0003 C 0004 C VARIABLES WHOSE DIMENSIONS ARE TO BE CHANGED 0005 C 0006 C RLAT (NSTA) JOIN(NSTA.,NJOINS) 0007 C RLONG (NSTA) P (MAXP ,MAXP ) 0008 C RSL (NSTA) T (NDEPTH,NWIDTH) 0009 C RCL (NSTA) A (MAXP ,MAXP2 ) 0010 C NOBS (NSTA) ITO (NJNM2) 0011 C IORDER(NSTA) ATA (NATA ) 0012 C ALONG(LASTRO) 0013 C 0014 C NUMBER(NSTA, 2) C (MAXP ) 0015 C NUMINV(LINV, 2) ICOL(MAXP2) 0016 C IASTAB(LASTRO,2) CON (MAXP2) 0017 C ELEV (NSTA, 2) DP (MAXP1) 0018 C NAME (NSTA, 4) X (NDEPTH) 0019 C TT (NSTA, 4) CK (NDEPTH) 0020 C 0021 C NSTA = TOTAL NO. OF STATIONS (FREE + FIXED) 0022 C MAXP = NUMBER OF CORRELATED OBSERVATIONS TO BE 0023 C HANDLED IN WEIGHT ARRAY P OR TWICE THE NO. 0024 C OF STATIONS FOR WHICH A POSITION WEIGHT 0025 C MATRIX IS REQUESTED 0026 C MAXP1= MAXP+1 0027 C MAXP2= 2*(MAXP+2) 0028 C NJOINS = MAX NUMBER OF CONNECTIONS TO ANY STATION 0029 C NJNM2 = NJOINS - 2 0030 C 0031 C IDEPTH = NUMBER OF FREE STATIONS 0032 C NDEPTH = IDEPTH*2 = DEPTH OF NORMALS 0033 C LWIDTH = SUM OF WIDTHS OF DIAG, VERT AND AUGMT 0034 C IN STATIONS 0035 C NWIDTH = LWIDTH * 2 = WIDTH OF NORMALS 0036 C LASTRD = MAX NO OF ASTRO CO-ORDINATES 0037 C NATA = MAXP2 * (MAXP2+1)/2 0038 C LINV = MAX NO OF INVERSE STATIONS 0039 C 0040 C ALL THE PARAMETERS FOR THE ABOVE VARIABLES MUST BE 0041 C CHOSEN AND/OR COMPUTED AND INSERTED IN THE 0042 C FOLLOWING DIMENSION STATEMENT 0043 C 0044 DIMENSION NUMBER( 99,2),NAME( 99,4),IORDER( 99), * RLONG( 99),RSL( 99),RCL( 99),NOBS( 99), * T(188,260),P(20,20),ALONG(15), * RLAT( 99),NUMINV(94,2),A(20,44), * ATA(990),C(20),ICOL(44), * DP(21),CK(188),X(188), * ELEV( 99,2),JOIN( 99,22),IASTAB(15,2), * CON(44),ITO(22),TT( 99,4) DOUBLE PRECISION RLAT, RLONG, ALONG, RSL, RCL, TT, C 0052 EQUIVALENCE (TT(1,1),T(1,1),JOIN(1,1)) C 0054 C THE VALUE OF THE FOLLOWING VARIABLE WORDS 0055 C MUST BE CHANGED SINCE THE ABOVE DIMENSION 0056 C STATEMENT WAS BASED ON THESE PARAMETERS 0057 C 0058 NSTA = 99 IDEPTH = 94 LWIDTH = 130 NJOINS = 22 LASTRO = 15 LINV = 94 MAXP = 20 C 0066 C THE FOLLOWING PARAMETERS ARE COMPUTED FROM THE ABOVE 0067 C BY THE PROGRAM 0068 C 0069 NDEPTH = IDEPTH * 2 0070 NWIDTH = LWIDTH * 2 0071 MAXP1 = MAXP + 1 0072 MAXP2 = (MAXP+2)*2 0073 NJNM2 = NJOINS - 2 0074 NATA = MAXP2*(MAXP2 + 1)/2 0075 CALL GAL1( NUMBER,NAME,IORDER,RLAT,RLONG,RSL, 0076 1 RCL, NOBS,JOIN,T,P,ALONG,IASTAB, 0077 2 NUMINV,A,ATA,C,CON,ICOL,DP,CK, 0078 3 X,ITO,TT,NSTA,IDEPTH,LWIDTH, 0079 4 NJOINS,MAXP,LASTRO,NDEPTH,NWIDTH, 0080 5 MAXP1,MAXP2,NJNM2,NATA,LINV,ELEV) 0081 STOP 0082 END 0083 SUBROUTINE INDEX (LINE, LPP, LINES, IHEAD ) 0086 DIMENSION IHEAD(20) 0087 IW=03 0131 IW=06 0131 IF(LINE-LPP+LINES) 10,10,20 0089 20 CONTINUE 0090 1018 WRITE (IW,1115) IHEAD 0091 1115 FORMAT (1H1,15X,40HGEODETIC SURVEY OF CANADA - PROGRAM GALS// 0092 1 10X, A3, 19A4 ) 0093 WRITE (IW,1019) 0094 1019 FORMAT (1H0,50X, 9HRESIDUALS/ 0095 1 1H0, 2X, 4HCODE, 2X, 12HSTATION FROM,12X, 10HSTATION TO, 0096 2 15X, 11HOBSERVATION, 6X, 3HWT , 5X, 3HRES, 6X, 3HPPM, 0097 1 4X, 9HADJD OBSN) 0098 LINE = 7 0099 RETURN 0100 10 LINE =LINE + LINES 0101 RETURN 0102 END 0103 SUBROUTINE GAL1( NUMBER,NAME,IORDER,RLAT,RLONG,RSL,RCL,NOBS,JOIN, 0106 1 T,P,ALONG,IASTAB,NUMINV,A,ATA,C,CON,ICOL,DP,CK, 0107 2 X,ITO,TT,NSTA,IDEPTH,LWIDTH,NJOINS,MAXP,LASTRO, 0108 3 NDEPTH,NWIDTH,MAXP1,MAXP2,NJNM2,NATA,LINV,ELEV) 0109 C GALS----GEOGRAPHICAL ADJUSTMENT BY LEAST SQUARES 0110 DIMENSION NUM(2), NAM(4), NUMFR(2),NUMTO(2), NAMFR(4), NAMTO(4), 0111 1 IHEAD(20), ICONTL(20) 0112 DIMENSION NUMBER(NSTA,2), NAME(NSTA,4), IORDER (NSTA), RLAT(NSTA), 0113 1 RLONG (NSTA), RSL(NSTA), RCL(NSTA), NOBS(NSTA), JOIN(NSTA, 0114 2 NJOINS), T(NDEPTH,NWIDTH), P(MAXP,MAXP), ALONG(LASTRO), 0115 3 IASTAB(LASTRO,2),NUMINV(LINV,2), A(MAXP,MAXP2), ATA(NATA), 0116 4 C(MAXP), CON(MAXP2), X (NDEPTH), DP(MAXP1), CK(NDEPTH), 0117 5 ICOL(MAXP2), ITO(NJNM2), TT(NSTA,4), ELEV(NSTA,2) 0118 DOUBLE PRECISION TT 0119 DOUBLE PRECISION RLAT,RLONG,ALONG,RSL,RCL 0120 DOUBLE PRECISIONDLAT,DMLAT,SLAT,DLONG,DMLONG, SLONG,DEG,DMIN,SECS, 0121 1OBS 0122 DOUBLE PRECISION DSIN,DCOS 0123 DOUBLE PRECISION C 0124 DATA NAM/4*4H / 0125 C WORDS IDEPTH,LASTRO,LWIDTH,NJOINS,MAXP 0126 C MUST BE CHANGED WHEN DIMENSIONS ARE CHANGED 0127 1 CONTINUE 0128 IR=05 0129 IW=06 0131 IP=07 IT1 = 8 0133 IT2 = 9 0134 IT3 =12 0135 REWIND IT1 0136 REWIND IT2 0137 C IT1,AND IT2 ARE TWO SCRATCH WORK AREAS 0138 C INPUT AND EDIT 0139 C USE INTEGER VARIABLE IPASS AS THE CARD COUNTER DURING EDITING 0140 5 FORMAT(1X,A3,19A4) 0141 15 FORMAT (1X, F9.5, I5, F5.2, 12I5) 0142 35 FORMAT(1X,I2,3X,6A4,1X,A2,5X,2(I4,I3,F9.5),2X,2A4) 0143 65 FORMAT ( 6H ERROR, I3, 22H - TOO MANY STATIONS (,I3, 1H)) 0144 85 FORMAT(1X,I2,3X,12A4,I4 ,I3 ,F9.5,F10.3) 0145 105 FORMAT(10F8.0) 0146 145 FORMAT ( 36H0CORRECT INPUT AS PER ERROR MESSAGES) 0147 165 FORMAT ( 6H0ERROR, I3, 25H - TOO MANY FREE STATIONS,2H (, I3, 1H)) 0148 185 FORMAT(6A4,I5,2G20.9) 0149 195 FORMAT (1H0///18X, 6HSTATNS/ 10X, 5HFIXED, I8/ 10X, 4HFREE, I9/ 0150 1 10X, 5HTOTAL,I8 / 10X, 9HDIAG BAND, I4/ 10X, 9HVERT BAND, I4/ 0151 2 10X, 8HAUG BAND, I5/ 10X, 5HASTRO, I8,/ 10X, 4HSPAN, I9) 0152 235 FORMAT(1H1,1X,A3,19A4) 0153 125 FORMAT(6H ERROR,I3,43H - INSUFFICIENT OBSERVATIONS TO FIX STATION, 0154 1 2X,2A4) 0155 C READ HEADING CARD 0156 READ(IR,5)IHEAD 0157 WRITE (IW,1115) IHEAD 0158 1115 FORMAT (1H1,15X,40HGEODETIC SURVEY OF CANADA - PROGRAM GALS// 0159 1 10X, A3, 19A4 ) 0160 WRITE (IW,1195) 0161 1195 FORMAT (1H0, 38X, 5HINPUT/) 0162 C READ CONTROL CARD 0163 READ (IR,15) CRITRN,NPASS,FICTSD, ISECTN, (ICONTL(I),I = 1,11) 0164 WRITE(IW,15) CRITRN,NPASS,FICTSD, ISECTN, (ICONTL(I),I = 1,11) 0165 NORMLS = ICONTL(1) 0166 C INPUT CO-ORDINATES 0167 C TEST FOR ERRORS IN PREVIOUS SECTION ADJUSTMENT 0168 IF(ISECTN) 2,3,2 0169 3 IERROR = 0 0170 2 IERROR = IERROR 0171 LSECT = IERROR 0172 IM = 1 0173 IA=0 0174 IASTRO=0 0175 IG=0 0176 IMCT=0 0177 IV=0 0178 N=0 0179 IPASS=0 0180 ISPAN=0 0181 IDIRFR = 0 0182 IWT = 0 0183 IFIX = 0 0184 C TEST FOR BACK SECTION ADJUSTMENT 0185 110 IF(ISECTN-2) 112,111,112 0186 111 IF(LSECT) 109,114,109 0187 114 READ(IT3) ICODE, NUM, SLAT, SLONG 0188 IF(ICODE-93) 113,109,113 0189 109 ISECTN = 0 0190 GO TO 112 0191 113 ICODE = 4 0192 IFIX = IFIX + 1 0193 WRITE (IW, 135) ICODE, NUM, SLAT, SLONG 0194 135 FORMAT ( 4X, I2, 3X, 2A4, 24X, 2G15.9 ) 0195 GO TO 72 0196 112 CONTINUE 0197 READ(IR,35)ICODE,NUM,NAM,IORD, LAT, MLAT,YLAT, LONG, MLONG,YLONG 0198 1 ,ELEV1,ELEV2 0199 WRITE(IW,35)ICODE,NUM,NAM,IORD, LAT, MLAT,YLAT, LONG, MLONG,YLONG 0200 1 ,ELEV1,ELEV2 0201 SLAT = YLAT 0202 SLONG = YLONG 0203 DMLONG = MLONG 0204 DLONG = LONG 0205 DMLAT = MLAT 0206 DLAT = LAT 0207 IPASS=IPASS+1 0208 IENTER = 3 0209 CALL INVALD (ICODE,IPASS,IERROR,IERR,IR,IW,IENTER) 0210 IF(IERR) 110,40,110 0211 40 IF(ICODE-10)30,20,10 0212 20 IFIX=IG 0213 GO TO 110 0214 30 CALL IGNOS (IENTER,IND,KEY,NUMBER,NUM ,NUM ,IG,IERROR, 0215 1 IR,IW, IPASS,IERR,IASTAB,IASTRO,ICODE,NSTA,LASTRO, 0216 1 LINV ) 0217 IF(IERR) 110,70,110 0218 70 CALL DMSR(DLAT,DMLAT,SLAT,SLAT,.5E-5) 0219 CALL DMSR(DLONG,DMLONG,SLONG,SLONG,.5E-5) 0220 IF(ICODE-7) 72,100,72 0221 72 IG=IG+1 0222 RLONG(IG) = SLONG 0223 RLAT(IG) = SLAT 0224 RSL(IG) = DSIN(RLAT(IG)) 0225 RCL(IG) = DCOS(RLAT(IG)) 0226 NUMBER(IG,1)=NUM(1) 0227 NUMBER(IG,2)=NUM(2) 0228 ELEV(IG,1) = ELEV1 0229 ELEV(IG,2) = ELEV2 0230 DO 71 J=1,4 0231 71 NAME(IG,J)=NAM(J) 0232 IORDER(IG)=IORD 0233 NOBS(IG)=0 0234 GO TO 110 0235 100 IASTRO=IASTRO+1 0236 IF(IASTRO-LASTRO) 101,101,102 0237 102 IERROR = 4 0238 WRITE (IW,55) IERROR,IASTRO 0239 55 FORMAT ( 6H ERROR,I3, 28H - TOO MANY ASTRO STATIONS (,I2,1H)) 0240 101 CONTINUE 0241 ALONG(IASTRO)= SLONG 0242 IASTAB(IASTRO,1) = NUM(1) 0243 IASTAB(IASTRO,2) = NUM(2) 0244 GO TO 110 0245 10 IF(ICODE-20)149,130,149 0246 130 IMCT=IG-IFIX 0247 GO TO 110 0248 149 IF(IMCT) 150,151,150 0249 151 IMCT = IG-IFIX 0250 150 IV=IG-IMCT-IFIX 0251 N=IG-IFIX 0252 140 IF(ICODE-30) 152,110,152 0253 152 IF(N-IDEPTH) 159,159,190 0254 190 IERROR = 10 0255 WRITE (IW,165) IERROR, N 0256 159 IF(NSTA-IG) 160,161,161 0257 160 IERROR = 11 0258 WRITE (IW, 65) IERROR, IG 0259 161 CONTINUE 0260 C END OF COORDINATE INPUT NOW START EDITING 0261 IPASS=0 0262 NFR=0 0263 NINV = 0 0264 DO 700 I=1,IG 0265 DO 700 J=1,NJOINS 0266 700 JOIN(I,J)=0 0267 600 READ(IR,85)ICODE,NUMFR,NAMFR,NUMTO,NAMTO,IEG, MIN,XECS,WT 0268 WRITE(IW,85)ICODE,NUMFR,NAMFR,NUMTO,NAMTO,IEG, MIN,XECS,WT 0269 SECS = XECS 0270 DEG = IEG 0271 DMIN = MIN 0272 200 CALL DMSR(DEG,DMIN,SECS,OBS,.5E-5) 0273 IPASS=IPASS+1 0274 IENTER = 2 0275 CALL INVALD (ICODE,IPASS,IERROR,IERR,IR,IW,IENTER) 0276 IF(IERR) 600,602,600 0277 602 IF(ICODE-80) 270,271,270 0278 270 IF(ICODE-81) 603,271,606 0279 606 IF(ICODE-82) 604,223,604 0280 271 CALL IABAND (IA,IPASS,IERROR,IT1,IW,IR,NUMBER,IG,IFIX, 0281 1 ICODE,NUMFR,IERR ,NOBS,IASTAB,IASTRO,NUMINV,NSTA,LASTRO, 0282 2 LINV) GO TO 600 0284 223 NINV = NINV + 1 0285 WRITE(IT2)ICODE,NUMFR,NAMFR,NUMTO,NAMTO,OBS,WT 0286 GO TO 600 0287 603 IF(ICODE-11) 605,604,605 0288 605 IENTER = 2 0289 CALL IGNOS (IENTER,IND,KEY,NUMBER,NUMFR,NUMTO,IG,IERROR, 0290 1 IR,IW, IPASS,IERR,IASTAB,IASTRO,ICODE,NSTA,LASTRO, 0291 2 LINV ) 0292 IF(IERR) 600,604,600 0293 604 CALL DIRSFR (ICODE,IW,NUMBER,NOBS, NAME ,IG,IDIRFR,KANG, 0294 1 KEY,IND,IFIX,IMCT,IM ,WT,IWT,IV,IA,IERROR,OBS,LWIDTH,MAXP, 0295 2 NSTA ) 0296 601 IF(ICODE-3)215,215,240 0297 215 IF(ICODE-2) 220,210,220 0298 210 OBS=SECS+DMIN*1000.+DEG*1000000. 0299 220 CONTINUE 0300 222 WRITE(IT1)ICODE,NUMFR,NAMFR,NUMTO,NAMTO,OBS,WT 0301 IF(ICODE-82) 250,600,250 0302 250 CALL LINK (JOIN, IND, KEY, NSTA, NJOINS, IERROR,NUMBER,IW) 0303 CALL LINK (JOIN, KEY, IND, NSTA, NJOINS, IERROR,NUMBER,IW) 0304 GO TO 600 0305 240 WRITE(IT1)ICODE,NUMFR,NAMFR,NUMTO,NAMTO,OBS,WT 0306 IF(ICODE-99)320,330,320 0307 320 IF(ICODE-98)340,330,340 0308 340 IF(ICODE-11) 351,360,351 0309 351 IF(LSECT)600,350,600 0310 350 CALL POSOBS(ICODE,IERROR,IW,IT1,IPASS,NOBS,NUMBER,IG,IR, 0311 1 IFIX,IMCT,IM,IV,IA,LWIDTH,MAXP,IASTAB,IASTRO,ISECTN,IT3,NSTA, 0312 2 P, LASTRO, LINV ) 0313 GO TO 600 0314 360 READ(IR,105) (( P(I,J),J=1,KANG),I = 1,KANG) 0315 WRITE(IW,105) (( P(I,J),J=1,KANG),I = 1,KANG) 0316 CALL TRIN ( P, KANG, 1, MAXP ) 0317 WRITE (IT1) (( P(I,J),J=1,KANG),I = 1,KANG) 0318 GO TO 600 0319 C CHECK IF SUFFICIENT OBSERVATIONS TO FIX STATION 0320 330 DO 500 J=1,IG 0321 IF(NOBS(J)-1)512,510,500 0322 510 IF(J-IFIX)500,500,512 0323 512 IERROR = 8 0324 WRITE (IW,125) IERROR,(NUMBER(J,K),K=1,2) 0325 500 CONTINUE 0326 WRITE (IW,195) IFIX, N, IG, IMCT,IV, IA, IASTRO, IM 0327 520 IF(IERROR)570,540,570 0328 570 WRITE(IW,145) 0329 IF(ICODE -98) 906,1,906 0330 540 CONTINUE 0331 IF(NINV) 544,544,542 0332 542 BACKSPACE IT1 0333 REWIND IT2 0334 REWIND IT3 0335 DO 543 I = 1,NINV 0336 READ (IT2)LCODE,NUMFR,NAMFR,NUMTO,NAMTO,OBS,WT 0337 543 WRITE(IT1)LCODE,NUMFR,NAMFR,NUMTO,NAMTO,OBS,WT 0338 WRITE(IT1)ICODE,NUMFR,NAMFR,NUMTO,NAMTO,OBS,WT 0339 544 CONTINUE 0340 REWIND IT1 0341 REWIND IT2 0342 REWIND IT3 0343 DO 900 I=1,IG 0344 900 WRITE(IT3)(JOIN(I,J),J=1,NJOINS) 0345 C DOUBLE UP NORMALS DIMENSION ON BASIS OF 2 UNKNOWNS PER STATION 0346 IM = 2*IM 0347 IV = 2*IV 0348 IA=2*IA 0349 N = 2*N 0350 C GAL2 SHOULD BE CALLED AT THIS POINT TO CARRY OUT THE ADJUSTMENT 0351 C THE VALUE OF CRIT SHOULD EVENTUALLY BE OPTIONALLY CONTROLLED 0352 CRIT = .4848E-9 0353 IF(CRITRN) 910,910,920 0354 920 CRIT = CRITRN*.485E-5 0355 910 CONTINUE 0356 CALL GAL2(RLAT,RLONG,RSL,RCL,ALONG,NUMBER,NAME,IASTAB,CRIT,IA,IG, 0357 1IFIX,IM,N,IASTRO,IV,IT1,IT2,IW,IR,IP,NPASS,MAXP,NUMINV,FICTSD,IT3, 0358 2 NORMLS, NSTA, NDEPTH, NWIDTH, LASTRO, MAXP1, MAXP2, NATA, 0359 3 A, ATA, C, CON, ICOL, DP, CK, X, T, P, IERROR, IHEAD, LINV ) 0360 L = 4 0361 REWIND IT3 0362 IF(ICONTL(2)) 904,905,904 0363 904 DO 907 I=1,IG 0364 READ (IT3) (NOBS(J), J = 1,NJOINS) 0365 907 TT(I,4) = I 0366 GO TO 903 0367 905 CONTINUE 0368 CALL GAL3(RLONG,RLAT,NUMBER,NAME,IORDER,IHEAD,IG,IW,IT3,RCL,RSL, 0369 1 NSTA, NDEPTH, NWIDTH, TT, T, ITO , NJNM2, ELEV ) 0370 903 CONTINUE 0371 WRITE (IW,1115) IHEAD 0372 WRITE (IW,1185) 0373 1185 FORMAT (1H0,14X, 40HINDEX AND ADJUSTED GEODETIC CO-ORDINATES/ 0374 1 1H0, 4HPAGE, 4X, 6HNUMBER, 5X, 4HNAME, 17X, 8HLATITUDE, 0375 2 9X, 9HLONGITUDE, 6X,9HELEVATION) 0376 DO 1171 J=1,IG 0377 CALL RDMS(DEG,DMIN,SECS,RLAT(J),.5E-5) 0378 IDEG = DEG 0379 IMIN = DMIN 0380 CALL RDMS(DEG,DMIN,SLAT ,RLONG(J),.5E-5) 0381 ID = DEG 0382 IIM= DMIN 0383 K = TT(J,4) 0384 WRITE (IW,1175) K,(NUMBER(J,K),K=1,2),(NAME(J,K),K=1,4), 0385 2 IDEG,IMIN, SECS, ID,IIM, SLAT , (ELEV(J,K),K=1,2) 0386 IF(ICONTL(3)) 908,909,908 0387 908 WRITE (IP, 35) L,(NUMBER(J,K),K=1,2),(NAME(J,K),K=1,4),IORDER(J), 0388 2 IDEG,IMIN, SECS, ID,IIM, SLAT , (ELEV(J,K),K=1,2) 0389 909 CONTINUE 0390 1175 FORMAT (1H ,I4,3X, 2A4, 3X, 4A4, 3X, 2I3, F8.4, 3X, 2I3, F8.4, 0391 1 5X, 2A4) 0392 1171 CONTINUE 0393 906 CONTINUE 0394 REWIND IT1 0395 REWIND IT2 0396 IF(ICODE -98) 571,1,571 0397 571 RETURN 0398 END 0399 SUBROUTINE INVALD (ICODE,IPASS,IERROR,IERR,IR,IW ,IENTER) 0402 C CHECKS VALIDITY OF CODES IN OBSERVATION DATA 0403 DIMENSION ICODES(20) 0404 DATA ICODES/1,2,3,80,92,81,91,93,98,99,11,82, 6,7,10,20,30,40,4,5/ 0405 IERR = 0 0406 K = 1 0407 L = 12 0408 IF(IENTER-3) 30,40,40 0409 40 K = 13 0410 L = 20 0411 30 DO 10 I = K,L 0412 IF(ICODE-ICODES(I)) 10,20,10 0413 10 CONTINUE 0414 IF(IENTER-3) 60,50,50 0415 50 IERR = 1 0416 GO TO 70 0417 60 IERR = 7 0418 70 IERROR = IERR 0419 WRITE (IW,45) IERR 0420 45 FORMAT ( 6H ERROR,I3,15H - INVALID CODE) 0421 20 RETURN 0422 END 0423 SUBROUTINE IGNOS (IENTER,IND,KEY,NUMBER,NUMFR,NUMTO,IG,IERROR, 0426 1 IR,IW, IPASS,IERR,IASTAB,IASTRO,ICODE,NSTA,LASTRO, 0427 2 LINV ) 0428 C CHECKS FOR PRESENCE OF CO-ORDINATES FOR OBSERVATIONS 0429 C ALSO CHECKS IF NUMBER FROM IS EQUAL TO NUMBER TO 0430 DIMENSION NUMBER(NSTA,2), NUMFR(2), NUMTO(2), IASTAB(LASTRO,2) 0431 IERR = 0 0432 IF(ICODE-7) 150,160,150 0433 160 CALL SEARCH (IND,IASTAB,NUMFR,IASTRO,LASTRO) 0434 GO TO 170 0435 150 CALL SEARCH ( IND, NUMBER, NUMFR, IG, NSTA ) 0436 170 IF(IND) 20,15,20 0437 15 IF(IENTER-2) 10,10,70 0438 10 IERR = 6 0439 IERROR = IERR 0440 55 FORMAT ( 6H ERROR,I3,29H - NO COORDINATE FOR STATION ,2A4) 0441 WRITE (IW,55) IERROR,(NUMFR ( J),J=1,2) 0442 20 IF(IENTER-2) 70,80,90 0443 90 IERR = 2 0444 IF(ICODE-7) 180,190,180 0445 190 IERR = 3 0446 180 IERROR = IERR 0447 WRITE (IW,45) IERROR 0448 45 FORMAT ( 6H ERROR,I3,22H - COORDINATE REPEATED) 0449 RETURN 0450 80 IF(ICODE-3) 100,110,100 0451 110 CALL SEARCH (INDAST,IASTAB,NUMFR,IASTRO,LASTRO ) 0452 IF (INDAST) 120,120,130 0453 120 IWARN = 4 0454 WRITE (IW,05) IWARN 0455 GO TO 100 0456 130 IWARN = 3 0457 WRITE (IW, 25) IWARN 0458 5 FORMAT ( 8H WARNING,I3, 27H - AZIMUTH ASSUMED GEODETIC) 0459 25 FORMAT ( 8H WARNING,I3, 29H - AZIMUTH ASSUMED ASTRONOMIC) 0460 100 CALL SEARCH (KEY,NUMBER,NUMTO,IG,NSTA) 0461 IF(KEY) 30,30,40 0462 30 IERR = 30 0463 IERROR = IERR 0464 WRITE (IW,55) IERROR,(NUMTO ( J),J=1,2) 0465 40 IF(KEY-IND) 70,50,70 0466 50 IERR = 5 0467 69 IERROR = IERR 0468 WRITE (IW,65) IERROR 0469 65 FORMAT ( 6H ERROR,I3,23H - SAME STATION NUMBERS) 0470 70 RETURN 0471 END 0472 SUBROUTINE IABAND (IA,IPASS,IERROR,IT1,IW,IR,NUMBER,IG,IFIX, 0475 1 ICODE,NUMFR,IERR ,NOBS,IASTAB,IASTRO,NUMINV,NSTA,LASTRO, 0476 2 LINV ) 0477 DIMENSION NUMBER(NSTA,2), NUMFR(2), NOBS(NSTA), IASTAB(LASTRO,2) 0478 1 ,NUMINV(LINV,2) 0479 45 FORMAT(8H WARNING,I3,42H - FIXED STATION IN AUGMENTED BAND IGNORED 0480 1 ) 0481 IF(ICODE - 80) 409,410,409 0482 410 IENTER = 1 0483 CALL IGNOS (IENTER,IND,KEY,NUMBER,NUMFR,NUMFR,IG,IERROR, 0484 1 IR,IW, IPASS,IERR,IASTAB,IASTRO,ICODE,NSTA,LASTRO, 0485 2 LINV ) 0486 IF(IERR) 429,414,429 0487 414 IF (IND - IFIX) 419,419,420 0488 419 IERR = 10 0489 C STATIONS IN THE AUGMENTED BAND CONNOT BE FIXED 0490 WRITE (IW,45) IERR 0491 RETURN 0492 420 NOBS (IND) = 1 0493 IA = IA + 1 0494 RETURN 0495 409 IF(IA) 430,429,430 0496 430 ICD80 = 80 0497 L = 1 0498 DO 20 I = 1,IG 0499 IF(NOBS (I)) 30,20,30 0500 30 WRITE (IT1) ICD80, (NUMBER(I,J),J = 1,2) 0501 NUMINV(L,1)=NUMBER(I,1) 0502 NUMINV(L,2)=NUMBER(I,2) 0503 L = L + 1 0504 20 NOBS (I) = 0 0505 ICD80 = 81 0506 WRITE (IT1) ICD80, (NUMBER(1,J),J = 1,2) 0507 IF(IA-LINV) 429,429,50 0508 50 IERROR = 26 0509 WRITE(IW,901) IERROR,IA 0510 901 FORMAT ( 6H ERROR, I3, 40H - TOO MANY INVERSE ELEMENTS REQUESTED ( 0511 1, I3, 1H)) 0512 429 RETURN 0513 END 0514 SUBROUTINE DIRSFR (ICODE,IW,NUMBER,NOBS, NAME ,IG,IDIRFR,KANG, 0517 1 KEY,IND,IFIX,IMCT,IM ,WT,IWT,IV,IA,IERROR,OBS,LWIDTH,MAXP, 0518 2 NSTA ) 0519 DIMENSION NUMBER(NSTA,2), NAME (NSTA,4), NOBS(NSTA) 0520 DOUBLE PRECISION TWOPI,OBS 0521 TWOPI = 3.1415926536D0 *2. 0522 MAXDIR = MAXP + 1 0523 IF(ICODE -1) 5,20,5 0524 5 NDIRFR = 0 0525 10 IF(IDIRFR-MAXDIR) 30,30,40 0526 40 IERROR = 15 0527 WRITE (IW,65) IERROR,IDIRFR,(NUMBER(INDTEM,J),J=1,2) 0528 65 FORMAT ( 6H ERROR,I3, 3H - ,21HTOO MANY DIRECTIONS (,I2, 0529 1 14H) FROM STATION,1X,2A4) 0530 GO TO 36 0531 20 IONE = 0 0532 IOBS = OBS/TWOPI 0533 OBS = OBS - IOBS*TWOPI 0534 IF(WT) 140, 130,140 0535 130 IONE = 1 0536 140 IF(IDIRFR) 50,50,60 0537 60 IWT = IWT + IONE 0538 IF(INDTEM-IND) 45,90,45 0539 90 IDIRFR = IDIRFR +1 0540 IENTER = 2 0541 C TEST FOR DIRECTION SEQUENCE 0542 ORIEND = OBS-ORIENT 0543 IF(ORIEND) 200,200,210 0544 200 ORIEND = ORIEND +TWOPI 0545 210 IF(ORIEND-OLDOBS) 220,220,230 0546 220 CONTINUE 0547 221 IERROR = 25 0548 WRITE (IW,225) IERROR 0549 225 FORMAT (6H ERROR,I3,3H - ,32HDIRECTION SEQUENCE OR REPETITION) 0550 GO TO 100 0551 230 OLDOBS = ORIEND 0552 100 CALL LSPAN (IENTER,IND ,KEY,IFIX,IMCT,IM,MIN,MAX,NOBS, 0553 1 IA,IV,IW,IERROR,LWIDTH, NSTA ) 0554 37 RETURN 0555 30 IF(IDIRFR - 1) 38,35,36 0556 38 IF(ICODE-3) 110,110,37 0557 C LENGTH AND AZIMUTH OBSERVATIONS WITH ZERO WEIGHTS WILL NOT BE 0558 C CONSIDERED IN THE SPAN COMPUTATION 0559 110 IF(WT) 39,120,39 0560 120 IERR = 21 0561 WRITE (IW,235) IERR 0562 235 FORMAT (8H WARNING,I3,34H - ZERO WEIGHT ON LAST OBSERVATION) 0563 RETURN 0564 35 IERROR = 16 0565 WRITE (IW,55) IERROR, (NUMBER(INDTEM,J),J=1,2) 0566 55 FORMAT(6H ERROR,I3,3H - ,32HONLY ONE DIRECTION FROM STATION ,2A4) 0567 36 KANG = IDIRFR - 1 0568 IF(IWT) 150,150,160 0569 160 IF(ICODE-11) 180,150,180 0570 180 IERROR = 18 0571 WRITE (IW,135) IERROR,IWT,(NUMBER(INDTEM,J),J=1,2) 0572 135 FORMAT ( 6H ERROR,I3,2H -,I3,1X, 0573 1 42HDIRECTIONS WITH ZERO WEIGHT OUT OF STATION,1X,2A4) 0574 150 IWT = 0 0575 IDIRFR = 0 0576 IF(NDIRFR) 50,38,50 0577 50 IDIRFR = IDIRFR +1 0578 IWT = IWT + IONE 0579 INDTEM = IND 0580 ORIENT = OBS 0581 KEYTEM = KEY 0582 OLDOBS = 0 0583 39 IENTER = 1 0584 GO TO 100 0585 45 NDIRFR = 1 0586 GO TO 10 0587 END 0588 SUBROUTINE LSPAN (IENTER,IND ,KEY,IFIX,IMCT,IM,MIN,MAX,NOBS, 0591 1 IA,IV,IW,IERROR,LWIDTH, NSTA ) 0592 C COMPUTES SPAN AND FILLS ARRAY NOBS 0593 DIMENSION NOBS(NSTA ) 0594 ISPAN = 0 0595 IF(IENTER-1) 300,300,620 0596 300 IF(IND-IFIX)301,301,610 0597 610 IF(IND-(IMCT+IFIX))621,621,301 0598 621 MAX=IND 0599 MIN=IND 0600 GO TO 620 0601 301 MAX = KEY 0602 MIN = KEY 0603 620 IF(KEY-IFIX)503,503,630 0604 630 IF(KEY-(IMCT+IFIX))302,302,503 0605 302 IF(KEY-MAX) 10,10,20 0606 10 IF(KEY-MIN)30,303,303 0607 30 MIN=KEY 0608 GO TO 40 0609 20 MAX=KEY 0610 IF(MIN) 40,30,40 0611 40 ISPAN=MAX-MIN+1 0612 IF(ISPAN-IM)303,303,315 0613 315 IM=ISPAN 0614 503 IF(IENTER-1) 303,403,303 0615 403 IF(KEY.GT.IFIX.AND.KEY.LE.(IMCT+IFIX)) GO TO 303 0616 603 MIN=0 0617 MAX=0 0618 303 NOBS(IND)=NOBS(IND)+1 0619 NOBS(KEY)=NOBS(KEY)+1 0620 IWIDTH = IA+IV+ISPAN 0621 IF(IWIDTH-LWIDTH) 50,50,60 0622 60 IERROR = 20 0623 WRITE (IW, 5) IERROR,ISPAN ,IWIDTH 0624 5 FORMAT(6H ERROR,I3,11H - TOO WIDE, 9H - SPAN =,I3, 9H, WIDTH =,I3) 0625 50 RETURN 0626 END 0627 SUBROUTINE POSOBS(ICODE,IERROR,IW,IT1,IPASS,NOBS,NUMBER,IG,IR, 0630 1 IFIX,IMCT,IM,IV,IA,LWIDTH,MAXP,IASTAB,IASTRO,ISECTN,IT3,NSTA, 0631 2 P, LASTRO, LINV ) 0632 DIMENSION NUM(2), NAM(4), NUMBER(NSTA,2), NOBS(NSTA), P(MAXP, 0633 1 MAXP), IASTAB(LASTRO,2) 0634 DOUBLE PRECISIONDLAT,DMLAT,SLAT,DLONG,DMLONG, SLONG,OBSLAT,OBSLG 0635 35 FORMAT(1X,I2,3X,6A4,A2,6X,2(I4 ,I3 ,F9.5)) 0636 45 FORMAT (6H ERROR,I3,16H - CODE SEQUENCE) 0637 105 FORMAT( 4E20.7) 0638 350 IF(ICODE-91)390,121,390 0639 121 IPE=0 0640 12 IF(ISECTN) 122,129,122 0641 122 READ (IT3)ICODE,NUM,OBSLAT,OBSLG 0642 WRITE (IW, 5) ICODE,NUM,OBSLAT,OBSLG 0643 5 FORMAT ( 1X, I2, 3X, 2A4, 24X, 2G15.9) 0644 GO TO 20 0645 129 READ(IR,35)ICODE,NUM,NAM,IORD, LAT, MLAT,XLAT, LONG, MLONG,XLONG 0646 WRITE(IW,35)ICODE,NUM,NAM,IORD, LAT, MLAT,XLAT, LONG, MLONG,XLONG 0647 SLONG = XLONG 0648 SLAT = XLAT 0649 DMLONG = MLONG 0650 DLONG = LONG 0651 DMLAT = MLAT 0652 DLAT = LAT 0653 CALL DMSR(DLAT,DMLAT,SLAT,OBSLAT,A) 0654 CALL DMSR(DLONG,DMLONG,SLONG,OBSLG,A) 0655 IENTER= 1 0656 CALL INVALD (ICODE,IPASS,IERROR,IERR,IR,IW,IENTER) 0657 IF(IERR) 12,20,12 0658 20 IF(ICODE-93)124,126,124 0659 126 WRITE(IT1)ICODE,NUM,OBSLAT,OBSLG 0660 IF(IPE) 390,390,127 0661 C NO POSITION OBSERVATIONS BETWEEN CODE 91 AND 93 0662 124 IF(ICODE - 92) 390,125,390 0663 C NO POSITION OBSERVATONS AFTER CODE 91 0664 125 IPE=IPE+2 0665 IENTER = 1 0666 CALL IGNOS (IENTER,IND,KEY,NUMBER,NUM ,NUM ,IG,IERROR, 0667 1 IR,IW, IPASS,IERR,IASTAB,IASTRO,ICODE,NSTA,LASTRO, 0668 2 LINV ) 0669 IF(IERR) 12,132,12 0670 132 NOBS(IND) = NOBS(IND) + 2 0671 WRITE(IT1)ICODE,NUM,OBSLAT,OBSLG 0672 IF(IPE-4) 400,410,420 0673 400 KEY = IND 0674 GO TO 12 0675 410 IENTER = 1 0676 GO TO 430 0677 420 IENTER = 3 0678 430 CALL LSPAN(IENTER,KEY,IND,IFIX,IMCT,IM,MIN,MAX,NOBS, 0679 1 IA,IV,IW,IERROR,LWIDTH, NSTA ) 0680 GO TO 12 0681 127 IF(IPE-MAXP) 198,198,210 0682 210 IERROR = 19 0683 IPEH = IPE/2 0684 WRITE(IW,55) IERROR,IPEH 0685 55 FORMAT ( 6H ERROR,I3, 31H - TOO MANY POSITION STATIONS (, I2, 1H)) 0686 198 IF(ISECTN) 199,200,199 0687 199 READ (IT3) ((P(I,J),J=1,IPE),I=1,IPE) 0688 GO TO 128 0689 200 READ(IR,105)((P(I,J),J=1,IPE),I=1,IPE) 0690 128 WRITE(IW,105)((P(I,J),J=1,IPE),I=1,IPE) 0691 WRITE(IT1) ((P(I,J),J=1,IPE),I=1,IPE) 0692 RETURN 0693 C ERROR IN CODE SEQUENCE 0694 390 IERROR=17 0695 WRITE (IW,45) IERROR 0696 RETURN 0697 END 0698 SUBROUTINE GAL3(RLONG,RLAT,NUMBER,NAME,IORDER,IHEAD,IG,IW,IT3,RCL, 0701 1RSL, NSTA, NDEPTH, NWIDTH, T, TT, ITO , NJNM2, ELEV ) 0702 DIMENSIONT(NSTA,4), TT(NDEPTH,NWIDTH), NUMBER (NSTA,2), NAME(NSTA, 0703 1 4), RLONG( NSTA), RLAT(NSTA), RSL (NSTA), RCL(NSTA), 0704 2 IORDER(NSTA), ITO(NJNM2), IHEAD(20), ELEV(NSTA,2) 0705 DOUBLE PRECISION RLONG,RLAT,RSL,RCL,ANGLE,TWOTAN,SA12,CA12,SA21, 0706 1CA21,R1,W1,R2,W2,S,AZD,AZM,AZS 0707 DOUBLE PRECISION T,TEMP 0708 C WE ARE USING 3 DOUBLE PRECISON COLUMS OF THE NORMALS ARRAY FOR SO 0709 REWIND IT3 0710 IPAGE=1 0711 DO 300 IC= 1, IG 0712 READ(IT3)IN,IFRM,(ITO(I),I=1,NJNM2) 0713 IFRM = IC 0714 DO 100 J=1,IN 0715 ITEMP=ITO(J) 0716 CALL GEOINV(RLONG(IFRM),RLONG(ITEMP),RSL(IFRM),RCL(IFRM),RSL(ITEMP 0717 1),RCL(ITEMP),R1,W1,R2,W2,SA12,CA12,SA21,CA21,S) 0718 ANGLE=TWOTAN(SA12,CA12) 0719 T(J,1) = ANGLE 0720 T(J,2) = S 0721 100 T(J,3) = ITEMP 0722 C SORT FOR INCREASING AZIMTUTHS AT THIS POINT 0723 IL = 1 0724 8 I=IN 0725 IF(I-IL)30,30,9 0726 9 J=I-1 0727 IF(I-IL)11,11,12 0728 11 IL = IL +1 0729 GO TO 8 0730 12 IF(T(I,1)-T(J,1))20,10,10 0731 10 I = I-1 0732 GO TO 9 0733 20 DO 25 K = 1,3 0734 TEMP = T(J,K) 0735 T(J,K) = T(I,K) 0736 25 T(I,K) = TEMP 0737 I = I-1 0738 GO TO 9 0739 30 CONTINUE 0740 C END OF SORT 0741 ICT=0 0742 LINECT=0 0743 DO 200 J = 1, IN 0744 ITEMP = T(J,3) 0745 IF(ICT) 69,69,50 0746 69 T(IC,4) = IPAGE 0747 GO TO 70 0748 50 IF(LINECT-18)80,80,60 0749 60 LINECT=0 0750 70 CALL TITLE(IHEAD,RLAT(IFRM),RLONG(IFRM),NUMBER(IFRM,1),NUMBER(IFRM 0751 1,2),IORDER(IFRM),IFRM,NAME,IW,IPAGE, ELEV(IFRM,1), ELEV(IFRM,2), 0752 1 NSTA ) 0753 80 CONTINUE 0754 ANGLE = T(J,1) 0755 S = T(J,2) 0756 CALL RDMS(AZD,AZM,AZS,ANGLE,5.E-03) 0757 IZD = AZD 0758 IZM=AZM 0759 IS = S + .0005 0760 Y=DABS(S-IS) WRITE(IW,85)IZD,IZM,AZS,IS,Y,(NAME(ITEMP,K),K=1,4), 0762 1 (NUMBER(ITEMP,K),K=1,2),IORDER(ITEMP) 0763 85 FORMAT(25X,2I4 ,F7.2,I12,F4.3,7X,3A4,A4,3X,A4,A4,4X,A2/) 0764 ICT=ICT+1 0765 LINECT=LINECT+1 0766 200 CONTINUE 0767 300 CONTINUE 0768 RETURN 0769 END 0770 SUBROUTINE GEOINV (WL1,WL2,SL1,CL1,SL2,CL2,R1,W1,R2,W2, 0773 1 SA12,CA12,SA21,CA21,S) 0774 DOUBLE PRECISION WL1,WL2,SL1,CL1,SL2,CL2,R1,W1,R2,W2,SA12,CA12, 0775 1SA21,CA21,S,A,E2,OME2,E2OME2,C1,C7,BR,SDL,CDL,SSL1,CA,OME2S2, 0776 2OME2S1,W2ON1,W2N1C2,CB,C1CA12,H2,BRH,C2CA21,WK,WKOR,F,F2,H, 0777 3C2,C3,C4,C5,C6,B ,BA 0778 DOUBLE PRECISION DSIN,DCOS,DSQRT 0779 C---- PARTIAL GEODETIC INVERSE 0780 C----- *INPUT, 0781 C---- WL1,WL2 = LONG. OF POINTS (1) AND (2) 0782 C---- SL1,CL1,SL2,CL2 = SINE,COSINE OF LAT. OF POINTS (1) AND (2) 0783 C----- *OUTPUT, 0784 C---- R1,W1,R2,W2 = RADII OF CURVATURE AT POINTS (1) AND (2) 0785 C---- SA12,CA12 = SINE,COSINE OF AZ. FROM (1) TO (2) 0786 C---- SA21,CA21 = SINE,COSINE OF AZ. FROM (2) TO (1) 0787 C---- S = COMPUTED DISTANCE FROM (1) TO (2) 0788 C----- 0789 A= 6378206.4 0790 B=6356583.8 0791 BA=B/A 0792 E2 = 1. - BA*BA 0793 OME2 = 1.-E2 0794 E2OME2 = E2/OME2 0795 C----- A= EQUATORIAL SEMI-AXIS 0796 C----- E2= EXCENTRICITY SQUARED 0797 C----- OME2= (1.-E2) 0798 C----- E2OME2= E2/(1.-E2) 0799 C----- 0800 C1 = .4166666667D-1 0801 C2=-.125 0802 C3= .0046875 0803 C4= .0375 0804 C5= .25 0805 C6= .1875 0806 C7 = .4166666667 D0 0807 BR= WL1-WL2 0808 SDL=DSIN (BR) 0809 CDL=DCOS (BR) 0810 SSL1= SL1*SL1 0811 W1= A/DSQRT(1.0-E2*SSL1) 0812 W2= A/DSQRT(1.0-E2*SL2*SL2) 0813 CA= CL2*SDL 0814 OME2S2= OME2*SL2 0815 CB= OME2S2*CL1-CL2*SL1*CDL+(W1/W2)*E2*CL1*SL1 0816 BR= 1.0/DSQRT(CA*CA+CB*CB) 0817 SA12= CA*BR 0818 CA12= CB*BR 0819 CA= -CL1*SDL 0820 OME2S1= OME2*SL1 0821 W2ON1= W2/W1 0822 W2N1C2= W2ON1*CL2 0823 CB= OME2S1*CL2-CL1*SL2*CDL+W2N1C2*E2*SL2 0824 BR= 1.0/DSQRT(CA*CA+CB*CB) 0825 SA21= CA*BR 0826 CA21= CB*BR 0827 C1CA12= CL1*CA12 0828 H2= E2OME2*C1CA12*C1CA12 0829 BRH= 1.0/(1.0+H2) 0830 R1= W1*BRH 0831 C2CA21= CL2*CA21 0832 R2= W2/(1.0+E2OME2*C2CA21*C2CA21) 0833 BR= W2N1C2*CDL-CL1 0834 CA= W2N1C2*SDL 0835 CB= W2ON1*OME2S2-OME2S1 0836 WK= W1*DSQRT(BR*BR+CA*CA+CB*CB) 0837 WKOR= WK/R1 0838 BR= E2OME2*SL1 0839 F= BR*C1CA12*BRH 0840 F2= F*F 0841 H= (BR*SL1-H2)*BRH 0842 S= WK*( 1. + WKOR*WKOR*( C1 + WKOR*( C2*F + WKOR*( C3 + C4*H 0843 1 + C5*F2 - WKOR*( C6*F*H + C7*F*F2 ) ) ) ) ) 0844 RETURN 0845 END 0846 SUBROUTINE SEARCH (IND,NUMBER,NUM,IG,NSTA) 0849 DIMENSION NUMBER(NSTA,2), NUM(2) 0850 IND=0 0851 IF(IG) 31,31,5 0852 5 DO 30 J=1,IG 0853 IF(NUMBER(J,1)-NUM(1))30,10,30 0854 10 IF(NUMBER(J,2)-NUM(2))30,20,30 0855 20 IND=J 0856 RETURN 0857 30 CONTINUE 0858 31 CONTINUE 0859 RETURN 0860 END 0861 SUBROUTINE DMSR(D,DM,S,R,ACC) 0864 DOUBLE PRECISION D,DM,S,R,RD,AM 0865 C CALL DMSR FOR DEGREES TO RADIANS 0866 C CALL RDMS FOR RADIANS TO DEGREES 0867 R = ((D*60.+DM)*60.+S)*.48481368111D-5 0868 GO TO 10 0869 ENTRY RDMS(D,DM,S,R,ACC) 0870 D=R*57.295779513D0 0871 KD=D 0872 RD = KD 0873 AM=(D-RD)*60. 0874 D=KD 0875 M=AM 0876 DM=M 0877 S=(AM-DM)*60. 0878 IF((S+ACC)-60.)2,1,1 0879 1 DM=DM+1 0880 S=0.0 0881 2 IF(DM-60.)10,3,3 0882 3 DM=0 0883 D=D+1 0884 10 CONTINUE 0885 RETURN 0886 END 0887 SUBROUTINE LINK (JOIN, IND, KEY, NSTA, NJOINS, IERROR,NUMBER,IW) 0890 DIMENSION JOIN(NSTA,NJOINS), NUMBER(NSTA,2) 0891 LIM=JOIN( IND,1)+2 0892 IF(LIM-2)70,70,50 0893 50 DO 40 J=3,LIM 0894 60 IF( JOIN(IND,J) - KEY) 40,20,40 0895 40 CONTINUE 0896 70 JOIN( IND,1)=JOIN( IND,1)+1 0897 JOIN(IND,LIM+1) = KEY 0898 IF(NJOINS-JOIN(IND,1))80,20,20 0899 80 IERROR = 35 0900 WRITE (IW,05) IERROR,JOIN(IND,1), (NUMBER(IND,I),I=1,2) 0901 5 FORMAT ( 6H ERROR, I3, 2H -, I4, 24H CONNECTIONS TO STATION ,2A4) 0902 20 RETURN 0903 END 0904 SUBROUTINE TITLE(IHEAD,PHI,AL,NUM,MUM ,IORD,INDENT,NAME,IW,I 0907 1PAGE, ELEV1, ELEV2, NSTA ) 0908 DIMENSION IHEAD(20), NAME(NSTA,4) 0909 DOUBLE PRECISION PHI,AL,PD,PM,PS,AD,AM,AS 0910 CALL RDMS(PD,PM,PS,PHI,5.E-05) 0911 CALL RDMS(AD,AM,AS,AL,5.E-05) 0912 WRITE(IW,5) IPAGE , IHEAD 0913 5 FORMAT (1H1,39X,40HGEODETIC SURVEY OF CANADA - PROGRAM GALS, 7X, 0914 1 4HPAGE, I5,///27X, A3, 19A4///) 0915 WRITE(IW,15) 0916 15 FORMAT(36X,7HSTATION,8X,4HNAME,4X,3HAND,5X,6HNUMBER, 0917 111X,5HORDER/) 0918 IPAGE = IPAGE + 1 0919 WRITE(IW,25)(NAME(INDENT,I),I=1,4),NUM,MUM,IORD 0920 25 FORMAT(48X,3A4,A4,3X,A4,A4, 10X, A2///) 0921 WRITE(IW,35) 0922 35 FORMAT (39X, 8HLATITUDE, 11X, 9HLONGITUDE, 9X, 9HELEVATION/) 0923 IPD = PD 0924 IPM = PM 0925 IAD=AD 0926 IAM=AM 0927 PY = PS 0928 AY = AS 0929 WRITE(IW,45)IPD,IPM,PY,IAD,IAM,AY, ELEV1, ELEV2 0930 45 FORMAT( 36X, 2(2I3, F8.4, 6X),2X,2A4//// 30X, 57HAZIMUTH AND D 0931 1ISTANCE TO STATION NAME AND NUMBER ,5X,5HORDER/) 0932 RETURN 0933 END 0934 DOUBLE PRECISION FUNCTION TWOTAN(A,B) 0937 C---- ARCTAN USING TWO ARGUMENTS 0938 C---- 0939 DOUBLE PRECISION A,B,TWOPI,PI,PIO2,AS,AC,W 0940 DOUBLE PRECISION DATAN 0941 TWOPI = 6.2831853072D0 0942 PI = 3.1415926536D0 0943 PIO2 = 1.5707963268D0 0944 AS =DABS (A) 0945 AC =DABS (B) 0946 IF(AS - AC) 20,20,30 0947 20 W = DATAN((AS/AC)) 0948 GO TO 40 0949 30 W =DATAN((AC/AS)) 0950 W =PIO2 - W 0951 40 IF(A) 50, 50, 80 0952 50 IF(B) 100, 100, 90 0953 80 IF(B) 70, 70 , 60 0954 70 W =PI - W 0955 GO TO 60 0956 90 W =TWOPI - W 0957 GO TO 60 0958 100 W =PI +W 0959 60 TWOTAN = W 0960 RETURN 0961 END 0962 SUBROUTINE GAL2(RLAT,RLONG,RSL,RCL,ALONG,NUMBER,NAME,IASTAB, 0965 1CRIT, IA,IG,IFIX,IM,N,IASTRO,IV,IWS,IRS,IW,IR,IP,NPASS,MAXP, 0966 2 NUMINV,FICTSD,IT3, NORMLS, NSTA, NDEPTH, NWIDTH, LASTRO, 0967 3 MAXP1, MAXP2, NATA, A, ATA,C, CON, ICOL, DP, CK, X, T, P, 0968 4 IERROR, IHEAD, LINV ) 0969 C GALS--GEOGRAPHICAL ADJUSTMENT BY LEAST SQUARES 0970 C ADJUSTMENT STAGE 0971 DIMENSION NUMBER(NSTA,2), NAME(NSTA,4), RLONG(NSTA), RLAT(NSTA), 0972 2 RSL(NSTA), RCL(NSTA), T(NDEPTH,NWIDTH), CK(NDEPTH), X(NDEPTH) 0973 2 , A(MAXP,MAXP2), P(MAXP,MAXP), ATA(NATA), C(MAXP), CON(MAXP2) 0974 3 , ICOL(MAXP2), IASTAB(LASTRO,2),ALONG(LASTRO), DP(MAXP1), 0975 4 AKI(5), STORE(5), ISFN(2), NAMFR(4), NUMTO(2), NAMTO(4), 0976 5 NUMINV(LINV,2), NUMFR(2) , IHEAD(20) 0977 DOUBLE PRECISION RLAT,RLONG,RSL,RCL,ALONG,C,AKI,STORE 0978 DOUBLE PRECISION ANGLE,OBS,R1,W1,R2,W2,SA12,CA12,SA21,CA21, 0979 1SCOMP,DUM,OBSLAT,OBSLG 0980 DOUBLE PRECISION DEG,DMIN,SECS 0981 DOUBLE PRECISION DSIN,DCOS 0982 IPASS=1 0984 VVSUM = 0 0985 F = 0. 0986 IRES = 0 0987 IFIRST = 0 0988 IF(NPASS) 1001,1001,1010 0989 1001 NPASS = 6 0990 1010 IT=IWS 0991 IWS=IRS 0992 IRS=IT 0993 LIM=IM+IV+IA 0994 DO 1011 I=1,N 0995 CK(I)=0.0 0996 DO 1011 J=1,LIM 0997 1011 T(I,J)=0.0 0998 1002 CONTINUE 0999 REWIND IWS 1000 REWIND IRS 1001 ICT = IM + IV 1002 IF(IA)1013,1013,1012 1003 C READ STATIONS TO BE AUGMENTED 1004 1012 READ(IRS)ICODE,NUMFR 1005 WRITE(IWS) ICODE,NUMFR 1006 IF(ICODE-80)1100,1014,1013 1007 C ASK FOR COV.TERMS ON FREE STATIONS ONLY,I.E. IF IFR IS LESS THAN 1008 C IFIX, WE HAVE AN ERROR CONDITION 1009 1014 IF(IRES)1009,1009,1008 1010 1008 WRITE(IW,5) ICODE,NUMFR 1011 GO TO 1012 1012 1009 CALL SEARCH(IFR,NUMBER,NUMFR,IG,NSTA) 1013 IF(IFR-IFIX)1100,1100,1015 1014 1015 ICT = ICT +1 1015 I = 2*(IFR-IFIX) -1 1016 T(I,ICT) = 1. 1017 ICT = ICT + 1 1018 T(I+1,ICT) = 1. 1019 GO TO 1012 1020 1013 NDIR=0 1021 KDIR=0 1022 IF(IRES)1020,1020,1018 1023 1018 LINE = 42 1024 CALL INDEX (LINE, 40, 7, IHEAD ) 1025 1020 READ(IRS)ICODE,NUMFR,NAMFR,NUMTO,NAMTO,OBS,WT 1026 WRITE(IWS)ICODE,NUMFR,NAMFR,NUMTO,NAMTO,OBS,WT 1027 IF(ICODE-1)1030,1200,1030 1028 1030 IF(ICODE-11)1040,1250,1040 1029 1040 IF(KDIR)1060,1060,1050 1030 1050 BACKSPACE IWS 1031 BACKSPACE IRS 1032 GO TO 1241 1033 1060 IF(ICODE-2)1070,1301,1070 1034 1070 IF(ICODE-3)1080,1300,1080 1035 1080 IF(ICODE-91)1090,1500,1090 1036 1090 IF(ICODE-99)1091,1110,1100 1037 C TO OBTAIN ERROR ELLIPSES 1038 1091 IF(ICODE-98)1101,1110,1100 1039 1101 IF ( ICODE - 82 ) 1100,1102,1100 1040 1102 IF(IRES) 1020,1020,1104 1041 1104 IF(IFIRST) 1109,1109,1103 1042 1103 CALL CURVE (IFIRST,T,NUMBER,NUMINV, RSL,RCL,X,NUMFR, 1043 1NUMTO,IA,IFIX,IG,ICODE,IV,IASTAB,IASTRO,IW,IPASS,IR,FICTSD,IM, 1044 2NSTA,NDEPTH, NWIDTH, IHEAD, LASTRO, LINV, RLONG, RLAT) 1045 GO TO 1020 1046 C END OF ERROR ELLIPSES 1047 C INPUT ERROR 1048 1100 WRITE(IW,45)ICODE 1049 45 FORMAT(28X,11HINPUT ERROR,I5) 1050 STOP 1051 1110 CONTINUE 1052 1109 IF(IFIRST) 9999,1131,9999 1053 1131 CONTINUE 1054 C IF IRES IS +VE, LIST S.E. OF UNIT WEIGHT 1055 IF(IRES) 1770,1770,1760 1056 1760 F= F-N 1057 SD = SQRT(VVSUM/F) 1058 SECS = SD/4.8481368111D-06 1059 SD = SD*1.E06 1060 XECS =SECS 1061 IF(FICTSD) 1762,1761,1762 1062 1761 IF(SD) 1764,1763,1764 1063 1763 FICTSD = 1. 1064 GO TO 1762 1065 1764 FICTSD = SD 1066 1762 CONTINUE 1067 IF = F 1068 WRITE (IW,1115) IHEAD 1069 1115 FORMAT (1H1,15X,40HGEODETIC SURVEY OF CANADA - PROGRAM GALS// 1070 1 10X, A3, 19A4 ) 1071 WRITE(IW,25) SD,XECS,IF 1072 25 FORMAT(1H0,10X,16HS.D. OF UNIT WT.,F10.4,2X,3HPPM,F10.4,2X,4HSECS/ 1073 1/10X,I10 ,2X,18HDEGREES OF FREEDOM) 1074 IF(ICODE-82) 9999,1103,9999 1075 C SUPRESS NORMALS VIA CONTROL CARD 1076 1770 IF(NORMLS) 1771,1334,1771 1077 1771 WRITE (IW,1115) IHEAD 1078 WRITE (IW,2015) 1079 2015 FORMAT (1H0, 50X, 7HNORMALS) 1080 WRITE (IW,1195) IPASS 1081 DO 1331 I=1,N 1082 WRITE(IW,1332)(T(I,J),J=1,LIM) 1083 1331 WRITE(IW,1333) CK(I) 1084 1332 FORMAT (1X, 6E10.3 ) 1085 1333 FORMAT (64X, E10.3 ) 1086 1334 CONTINUE 1087 WRITE (IW,1115) IHEAD 1088 WRITE (IW,1195) IPASS 1089 1195 FORMAT (5H0PASS,I3/) 1090 CALL CHOLES(T,CK,X,IM,IV,IA,N, NDEPTH, NWIDTH ) 1091 RS = 1./.48481368111D-5 1092 ABST=0.0 1093 J = IFIX 1094 IF(NPASS-1) 1120,1120,1119 1095 1119 CONTINUE 1096 WRITE (IW,2175) 1097 2175 FORMAT ( 5H0PAGE,4X,6HNUMBER, 5X, 4HNAME, 1098 2 17X, 8HLATITUDE,9X, 5HCORRN, 8X, 9HLONGITUDE,8X, 5HCORRN) 1099 DO 1113 I = 1, N,2 1100 J = J+1 1101 IF(ABS(X(I+1))-ABS(X(I))) 1116,1116,1117 1102 1117 CHK = X(I+1) 1103 GO TO 1118 1104 1116 CHK = X(I) 1105 1118 IF(ABS(CHK)-ABST) 1112,1112,1111 1106 1111 ABST = ABS(CHK) 1107 1112 RLAT(J)=RLAT(J) + X(I) 1108 RLONG(J)=RLONG(J)+X(I+1) 1109 CALL RDMS(DEG,DMIN,SECS,RLAT(J),.5E-5) 1110 IDEG = DEG 1111 IMIN = DMIN 1112 CALL RDMS(DEG,DMIN,SCOMP,RLONG(J),.5E-5) 1113 ID = DEG 1114 IIM= DMIN 1115 DEG = X(I)*RS 1116 DMIN = X(I+1)*RS 1117 WRITE (IW,55) J, (NUMBER(J,K),K=1,2),(NAME(J,K),K=1,4), 1118 1 IDEG,IMIN, SECS, DEG, ID,IIM, SCOMP, DMIN 1119 55 FORMAT (1H ,I4,3X, 2A4, 3X, 4A4, 3X, 2I3, F8.4, 2X, F10.4, 4X, 1120 1 2I3, F8.4, 2X, F10.4) 1121 1113 CONTINUE 1122 DO 1114 L=1,IG 1123 RSL(L) =DSIN(RLAT(L)) 1124 RCL(L) =DCOS (RLAT(L)) 1125 1114 CONTINUE 1126 C CRIT IS A VALUE FROM THE CONTROL CARD 1127 IF(IPASS-2) 1120,1138,1138 1128 1138 IF(CRIT-ABST) 1120,1120,1140 1129 1120 IPASS=IPASS+1 1130 IF(IPASS-NPASS) 1010,1010,1121 1131 C CONVERGENCE HAS NOT BEEN REACHED AFTER 3 PASSES 1132 1121 WRITE(IW,1122) 1133 IERROR = 40 1134 1122 FORMAT(10X,26HSOLUTION HAS NOT CONVERGED) 1135 IF(FICTSD) 1139,1139,1140 1136 1139 FICTSD = 1. 1137 GO TO 1140 1138 1200 IF(KDIR)1210,1210,1220 1139 1210 KDIR=1 1140 NDIR=0 1141 ISFN(1)=NUMFR(1) 1142 ISFN(2)=NUMFR(2) 1143 CALL SEARCH(IFR,NUMBER,NUMFR,IG,NSTA) 1144 IF(IFR-IFIX) 1212,1212,1213 1145 1212 ICOL(1) = 0 1146 ICOL(2) = 0 1147 GO TO 1214 1148 1213 ICOL(1) = 2*(IFR-IFIX) - 1 1149 ICOL(2)= ICOL(1)+ 1 1150 1214 JCOL = 0 1151 IROW =0 1152 DO 1211 I=1,MAXP 1153 DO 1211 J=1,MAXP2 1154 1211 A(I,J)=0.0 1155 GO TO 1230 1156 1220 IF(NUMFR(2)-ISFN(2))1050,1221,1050 1157 1221 IF(NUMFR(1)-ISFN(1))1050,1230,1050 1158 1230 CALL SEARCH(ITO,NUMBER,NUMTO,IG,NSTA) 1159 C IF IRES IS POSITIVE LIST DIRECTIONS 1160 IF(IRES)1880,1880,1870 1161 1870 CALL RDMS(DEG,DMIN,SECS,OBS,5.E-04) 1162 CALL INDEX (LINE,40,1,IHEAD) 1163 ID= DEG 1164 MIN = DMIN 1165 XECS = SECS 1166 WRITE(IW,5)ICODE,NUMFR,NAMFR,NUMTO,NAMTO,ID,MIN,XECS,WT 1167 5 FORMAT(1X,I5,12A4,I5 ,I3,F7.3,2F 9.3) 1168 1880 IF(ITO-IFIX)1215,1215,1216 1169 1215 ICOL(JCOL+3) = 0 1170 ICOL(JCOL+4) = 0 1171 GO TO 1217 1172 1216 ICOL(JCOL+3) = 2*(ITO-IFIX) -1 1173 ICOL(JCOL+4)=ICOL(JCOL+3) +1 1174 1217 NDIR = NDIR + 1 1175 C ARBITRARILY WE SET ALL ZERO WTS. = 0.001 1176 IF(WT) 1219,1218,1219 1177 1218 WT=0.001 1178 1219 DP(NDIR)=WT 1179 L=1 1180 CALL GEOINV(RLONG(IFR),RLONG(ITO),RSL(IFR),RCL(IFR),RSL(ITO), 1181 1RCL(ITO),R1,W1,R2,W2,SA12,CA12,SA21,CA21,SCOMP) 1182 CALL OBSEQN (RCL(IFR),RCL(ITO),R1,W1,R2,W2,SA12,CA12,SA21,CA21, 1183 1OBS,SCOMP,1,L,AKI,RSL(IFR),RLONG(IFR),DUM) 1184 IF(NDIR-1)1233,1233,1232 1185 1232 A(IROW,1) =AKI(1)-STORE(1) 1186 A(IROW,2) =AKI(2)-STORE(2) 1187 A(IROW,JCOL+1)= -STORE(3) 1188 A(IROW,JCOL+2) = -STORE(4) 1189 A(IROW,JCOL+3) = AKI(3) 1190 A(IROW,JCOL+4) = AKI(4) 1191 C(IROW) = AKI(5)- STORE(5) 1192 C(IROW) = ANGLE(C(IROW)) 1193 1233 DO 1231 I=1,5 1194 1231 STORE(I)=AKI(I) 1195 JCOL=JCOL+2 1196 IROW=IROW+1 1197 GO TO 1020 1198 1241 KANG= NDIR -1 1199 DO 1243 I=1,KANG 1200 DO 1242 J=1,KANG 1201 1242 P(I,J)=0.0 1202 1243 P(I,I)=1./DP(I) 1203 DO 1245 I=1,KANG 1204 P(I,I)=P(I,I)+1./DP(I+1) 1205 IF(I+1-KANG) 1244,1244,1245 1206 1244 P(I,I+1)= -1./DP(I+1) 1207 P(I+1,I)= -1./DP(I+1) 1208 1245 CONTINUE 1209 CALL TRIN (P, KANG, 1, MAXP ) 1210 ICODE=11 1211 WRITE(IWS)ICODE ,NUMFR,NAMFR,NUMTO,NAMTO,OBS,WT 1212 WRITE(IWS)((P(I,J),J=1,KANG),I=1,KANG) 1213 GO TO 999 1214 1250 KANG =NDIR-1 1215 READ(IRS)((P(I,J),J=1,KANG),I=1,KANG) 1216 WRITE(IWS)((P(I,J),J=1,KANG),I=1,KANG) 1217 C FORM NORMALS CONTRIBUTION NOW 1218 999 IROW=IROW-1 1219 JCOL = JCOL + 2 1220 C WHEN IRES IS POSITIVE , LIST RESIDUALS 1221 IF(IRES)1998,1998,1997 1222 1997 DO 1993 I = 1, KANG 1223 1993 CK(I)=-C(I)/4.8481368111D-06 1224 LINES = KANG/8 + 4 1225 CALL INDEX ( LINE,40,LINES,IHEAD) 1226 WRITE(IW,35) (CK(I),I=1,KANG) 1227 35 FORMAT( 5X,15HANGLE RESIDUALS/5X,( 8F10.3)) 1228 WRITE(IW,36) 1229 36 FORMAT(/) 1230 1998 CALL MATMUL(P,A,C,KANG,JCOL, MAXP, MAXP2) 1231 C OBSERVATIONS ARE WEIGHTED NOW 1232 IF(IRES) 1996,1996,1995 1233 1995 DO 1994 I = 1, KANG 1234 1994 VVSUM = VVSUM + C(I)*C(I) 1235 F = F + KANG 1236 KDIR = 0 1237 NDIR = 0 1238 GO TO 1020 1239 1996 DO 1640 I = 1, JCOL 1240 DO 1620 K=I,JCOL 1241 SUM=0.0 1242 DO 1610 J=1,IROW 1243 1610 SUM=SUM+ A(J,I)*A(J,K) 1244 IK = (I-1)*JCOL -(I*I-I)/2 + K 1245 1620 ATA(IK)=SUM 1246 SUM=0.0 1247 DO 1630 J=1,IROW 1248 1630 SUM =SUM +A(J,I)*C(J) 1249 1640 CON(I)=SUM 1250 1330 CALL TRAN (ATA,CON,ICOL,T,CK,JCOL,IM,IV,IA,N,NDEPTH, NWIDTH, NATA, 1251 1 MAXP2) 1252 KDIR=0 1253 NDIR=0 1254 GO TO 1020 1255 1300 CALL SEARCH(IND,IASTAB,NUMFR,IASTRO, LASTRO ) 1256 IF(IND)1301,1301,1308 1257 1308 DUM=ALONG(IND) 1258 1301 CALL SEARCH(IFR,NUMBER,NUMFR,IG, NSTA) 1259 IF(WT) 1312,1020,1312 1260 1312 CONTINUE 1261 CALL SEARCH(ITO,NUMBER,NUMTO,IG,NSTA) 1262 IF(IND) 1309,1309,1311 1263 1309 DUM = RLONG(IFR) 1264 1311 IF(IFR-IFIX)1302,1302,1303 1265 1302 ICOL(1) =0 1266 ICOL(2)=0 1267 GO TO 1306 1268 1303 ICOL(1) = 2*(IFR-IFIX) -1 1269 ICOL(2)=ICOL(1)+1 1270 1306 IF(ITO-IFIX)1304,1304,1305 1271 1304 ICOL(3) =0 1272 ICOL(4)=0 1273 GO TO 1307 1274 1305 ICOL(3)=2*(ITO-IFIX) - 1 1275 ICOL(4)=ICOL(3)+1 1276 1307 L=1 1277 CALL GEOINV(RLONG(IFR),RLONG(ITO),RSL(IFR),RCL(IFR),RSL(ITO), 1278 1RCL(ITO),R1,W1,R2,W2,SA12,CA12,SA21,CA21,SCOMP) 1279 CALL OBSEQN(RCL(IFR),RCL(ITO),R1,W1,R2,W2,SA12,CA12,SA21,CA21, 1280 1OBS,SCOMP,ICODE,L,AKI, RSL(IFR),RLONG(IFR),DUM) 1281 IF(ICODE-3) 1319,1318,1319 1282 1318 CONTINUE 1283 1319 CONTINUE 1284 AKI(5) = ANGLE(AKI(5)) 1285 C IF IRES IS POSITIVE, LIST RESIDUALS AND CALCULATE S.D. 1286 IF(IRES) 1790,1790,1730 1287 1730 VVSUM = VVSUM + AKI(5)*WT*AKI(5) 1288 CALL INDEX (LINE,42, 1, IHEAD ) 1289 F = F+1 1290 IF(ICODE-2)1750,1720,1750 1291 1720 PPM =-AKI(5) * 1.E06 1292 AKI(5) = -SCOMP*AKI(5) 1293 IOBS = OBS + .0005 1294 XBS=DABS(OBS-IOBS) ICOMP = SCOMP + .0005 1296 XCOMP=DABS(SCOMP-ICOMP) WRITE(IW,1735)ICODE,NUMFR,NAMFR,NUMTO,NAMTO, 1298 1 IOBS, XBS, WT, AKI(5), PPM, ICOMP, XCOMP 1299 1735 FORMAT (1X, I5, 12A4, I11, F4.3, 3F9.3, I9, F4.3) 1300 GO TO 1020 1301 1750 CALL RDMS(DEG,DMIN,SECS,OBS,5.E-04) 1302 ID = DEG 1303 MIN = DMIN 1304 AKI(5) = -AKI(5)/4.8481368111D-06 1305 WRITE(IW,5) ICODE,NUMFR,NAMFR,NUMTO,NAMTO,ID, MIN,SECS,WT,AKI(5) 1306 GO TO 1020 1307 1790 CONTINUE 1308 DO 1320 I=1,4 1309 DO 1310 K=I,4 1310 IK= (I-1)*4 -(I*I-I)/2 +K 1311 1310 ATA(IK) =AKI(I)*WT*AKI(K) 1312 1320 CON(I)=AKI(I)*WT*AKI(5) 1313 JCOL=4 1314 GO TO 1330 1315 C POSITION EQUATION INPUT 1316 1500 ICT =0 1317 IOBS=1 1318 1501 READ(IRS)ICODE,NUMFR,OBSLAT,OBSLG 1319 WRITE(IWS)ICODE,NUMFR,OBSLAT,OBSLG 1320 IF(ICODE-92)1503,1502,1503 1321 1502 CALL SEARCH(IFR,NUMBER,NUMFR,IG,NSTA) 1322 ICT=IOBS 1323 IF(IFR-IFIX)1504,1504,1505 1324 1504 ICOL(ICT) =0 1325 ICOL(ICT+1) =0 1326 GO TO 1506 1327 1505 ICOL(ICT)=2*(IFR-IFIX) - 1 1328 ICOL(ICT+1) = ICOL(ICT)+1 1329 1506 C(IOBS)=OBSLAT - RLAT(IFR) 1330 C(IOBS+1)= OBSLG-RLONG(IFR) 1331 IOBS = IOBS + 2 1332 GO TO 1501 1333 1503 IOBS=IOBS -1 1334 READ(IRS)((P(I,J),J=1,IOBS), I=1,IOBS) 1335 WRITE(IWS)((P(I,J),J=1,IOBS),I=1,IOBS) 1336 DO 1520 I=1,IOBS 1337 SUM=0 1338 DO 1510 K=I,IOBS 1339 IK = (I-1)*IOBS -(I*I-I)/2 + K 1340 ATA(IK) =P(I,K) 1341 1510 CONTINUE 1342 DO 1511 L = 1,IOBS 1343 1511 SUM=SUM + P(I,L)*C(L) 1344 1520 CON(I) = SUM 1345 JCOL =IOBS 1346 GO TO 1330 1347 C RESIDUAL AND LISTING STAGE 1348 1140 CONTINUE 1349 IRES = 1 1350 C PUT OUT PERTINENT PARTS OF THE INVERSE AT THIS POINT 1351 IF(IA) 1010,1010,1144 1352 1144 LLIM = LIM - IA + 1 1353 WRITE(IW,1145) 1354 1145 FORMAT(1H1,10X,16HINVERSE ELEMENTS/) 1355 DO 1143 I = 1, IA 1356 WRITE(IW,1146) (T(I,J),J=LLIM,LIM) 1357 1143 LLIM =LLIM + 1 1358 1146 FORMAT( 4E20.7) 1359 C OBTAIN WEIGHT MATRIX BY INVERTING AUGMENTED BAND ELEMENTS 1360 IF(IA-MAXP) 1149,1149,1148 1361 1148 WRITE (IW,1155) 1362 1155 FORMAT (46H0AUGMENTED BAND MATRIX TOO LARGE FOR INVERSION) 1363 IERROR = 41 1364 GO TO 1002 1365 1149 WRITE (IW,1165) 1366 1165 FORMAT(////10X,22HPOSITION WEIGHT MATRIX/) 1367 J = LIM-IA+1 1368 DO 1150 I=1,IA 1369 DO 1147 JN=I,IA 1370 P(I,JN) = T(I,J) 1371 P(JN,I) = T(I,J) 1372 1147 J = J+1 1373 1150 J = LIM-IA+1+I 1374 CALL TRIN(P,IA,2,MAXP) 1375 K = IA/2 1376 IK = 91 1377 WRITE (IW,905) IK 1378 905 FORMAT ( 1X,I2 ) 1379 IK = 92 1380 DO 908 I = 1,K 1381 NUMTO(1) = NUMINV(I,1) 1382 NUMTO(2) = NUMINV(I,2) 1383 CALL SEARCH (IND, NUMBER, NUMTO, IG, NSTA ) 1384 WRITE (IW,915) IK,(NUMBER(IND,J),J=1,2), RLAT(IND),RLONG(IND) 1385 915 FORMAT ( 1X, I2, 3X, 2A4, 24X, 2G15.9) 1386 908 WRITE (IT3) IK,(NUMBER(IND,J),J=1,2) , RLAT(IND), RLONG(IND) 1387 IK = 93 1388 WRITE (IT3) IK,(NUMBER(IND,J),J=1,2) , RLAT(IND), RLONG(IND) 1389 WRITE (IW,905) IK 1390 WRITE (IT3) ((P(I,J),J=1,IA), I=1,IA) 1391 DO 1151 I = 1, IA 1392 1151 WRITE(IW,1146)(P(I,J) , J = 1, IA) 1393 WRITE (IP,1146)((P(I,J),J=1,IA),I=1,IA) 1394 GO TO 1002 1395 9999 CONTINUE 1396 RETURN 1397 END 1398 SUBROUTINE TRAN (ATA,CON,ICOL,T,C,JCOL,IM,IV,IA,N, NDEPTH, NWIDTH, 1401 1 NATA, MAXP2 ) 1402 DIMENSION ATA(NATA), T(NDEPTH,NWIDTH), C(NDEPTH), CON(MAXP2), 1403 1 ICOL(MAXP2) 1404 C THIS SUBROUTINE REQUIRES FUNCTION JSUB ASSOCIATED WITH CHOLES 1405 DO 20 I=1,JCOL 1406 IF(ICOL(I))20,20,7 1407 7 ISUB= ICOL(I) 1408 C(ISUB) = C(ISUB) + CON(I) 1409 DO 10 J = I,JCOL 1410 IF(ICOL(J)) 10,10,9 1411 9 IJ=(I-1)*JCOL -(I*I-I)/2 + J 1412 IC = ICOL(I) 1413 JCT = ICOL(J) 1414 IF(IC-JCT)11,11,8 1415 8 IC = ICOL(J) 1416 JCT = ICOL(I) 1417 11 JC = JSUB(IC,JCT,IM,IV,IA,N) 1418 T(IC,JC) = ATA(IJ)+T(IC,JC) 1419 10 CONTINUE 1420 20 CONTINUE 1421 RETURN 1422 END 1423 SUBROUTINE TRIN(A,M,IENTER,MAXP) 1426 DIMENSION A(MAXP,MAXP) 1427 IF(M-1)3,3,200 1428 3 A(1,1) = SQRT(1./A(1,1)) 1429 GO TO 300 1430 200 DO 7 L=1,M 1431 FMN=1./A(L,1) 1432 DIV=A(L,1) 1433 DO 4 J=2,M 1434 4 A(L,J-1)=A(L,J)/DIV 1435 A(L,M) =FMN 1436 DO 7 J=1,M 1437 IF(J-L)5,7,5 1438 5 FMULT=A(J,1) 1439 DO 6 K=2,M 1440 6 A(J,K-1)= A(J,K)-FMULT*A(L,K-1) 1441 A(J,M)= -FMN*FMULT 1442 7 CONTINUE 1443 IF(IENTER-1) 8,8,300 1444 8 CONTINUE 1445 A(1,1) = SQRT (A(1,1)) 1446 DO 10 I = 2,M 1447 10 A(1,I) = A(1,I)/A(1,1) 1448 DO 100 I = 2,M 1449 DIV = 0. 1450 N = I-1 1451 DO 30 L = 1, N 1452 30 DIV = DIV + A(L,I)*A(L,I) 1453 A(I,I) = SQRT (A(I,I) - DIV) 1454 JL = I + 1 1455 IF(JL-M) 20,20,100 1456 20 DO 90 J = JL,M 1457 DIV = 0. 1458 DO 40 L = 1,N 1459 40 DIV = A(L,I)*A(L,J) + DIV 1460 A(I,J) = (A(I,J) - DIV)/A(I,I) 1461 90 CONTINUE 1462 100 CONTINUE 1463 300 CONTINUE 1464 RETURN 1465 END 1466 SUBROUTINE MATMUL(P,A,C,N,M, MAXP, MAXP2 ) 1469 DIMENSION P(MAXP,MAXP), A(MAXP,MAXP2), C(MAXP) 1470 C MATRIX C IS THE CONSTANTS AND IS N BY 1 1471 C MATRIX P IS UPPER TRIANGULAR AND IS N BY N 1472 C MATRIX A IS N BY M 1473 DOUBLE PRECISION C 1474 DO 20 I=1,N 1475 DO 20 K=1,M 1476 SUM=0.0 1477 DO 10 J=I,N 1478 10 SUM=SUM +P(I,J)*A(J,K) 1479 20 A(I,K)=SUM 1480 DO 40 I = 1, N 1481 SUM = 0. 1482 DO 30 J = I,N 1483 30 SUM = SUM +P(I,J)*C(J) 1484 40 C(I) = SUM 1485 RETURN 1486 END 1487 DOUBLE PRECISION FUNCTION ANGLE(CONST) 1490 DOUBLE PRECISION CONST 1491 IF(DABS(CONST)-1.570796327D0)4,4,1 1492 1 IF(CONST) 2,2,3 1493 2 CONST = CONST +6.283185307D0 1494 GO TO 4 1495 3 CONST = CONST -6.283185307D0 1496 4 ANGLE = CONST 1497 RETURN 1498 END 1499 SUBROUTINE OBSEQN (CL1,CL2,R1,W1,R2,W2,SA12,CA12,SA21,CA21,OBS,S, 1502 1 K, L , AKI, SL1,WL,ASL) 1503 DOUBLE PRECISION AKI,CL1,CL2,R1,W1,R2,W2,SA12,CA12,SA21,CA21,OBS, 1504 1S,SL1,WL,ASL,SLA,OS,APR,TWOTAN 1505 C---- FORM OBSERVATION EQUATION 1506 C** INPUT 1507 C SL1 = SINE OF LAT. OF POINT 1 1508 C WL = APPROX. LONG. OF POINT 1 1509 C ASL = ASTRO. LONG. OF POINT 1 1510 C---- CL1, CL2 =COSINE OF LAT. OF POINTS 1 AND 2 1511 C SA12 , CA12 =SIN AND COS OF AZ FROM 1 TO 2 1512 C SA21 , CA21 = ,, ,, A2 ,, 2 ,, 1 1513 C OBS = OBSERVATION , S= COMP. LENGTH 1514 C K = 2 ,, IS LENGTH 1515 C---- K = 3 ,, IS AZ 1516 C---- K = 1 ,, IS DIRECTION 1517 C L = 2 COMPUTE RESIDUALS 1518 C L = 1 NORMAL 1519 C** OUTPUT 1520 C---- AKI(1 TO 4) = FOUR ELEMENTS OF OBS. EON. 1521 C---- AKI(5) = CONSTANT TERM, WHEN L=2, AKI(5) = RESIDUAL 1522 DIMENSION AKI( 5) 1523 SLA = SL1 1524 GO TO (300,100,200),K 1525 100 GO TO (101, 102),L 1526 101 OS = 1.0 / S 1527 AKI(1) = OS *(-R1) * CA12 1528 AKI(2) = OS * W1 * CL1 * SA12 1529 AKI(3) = OS *(-R2) * CA21 1530 AKI(4) = OS * W2 * CL2 * SA21 1531 102 AKI(5)=OBS*OS-1. 1532 RETURN 1533 300 SLA = 0. 1534 200 GO TO (201, 202),L 1535 201 OS = 1. / S 1536 AKI(1) = OS * SA12* R1 1537 AKI(2) = OS * CA12 * CL1*W1 AKI(3) = OS * SA21*R2 1539 AKI(4) = OS * CA21 * CL2* W2 1540 202 APR = TWOTAN(SA12, CA12) 1541 AKI(5) = OBS - APR -(WL - ASL)*SLA 1542 RETURN 1543 END 1544 SUBROUTINE CHOLES(T,C,X,IM,IV,IA,N,NDEPTH, NWIDTH ) 1547 DIMENSION T(NDEPTH,NWIDTH), C(NDEPTH), X(NDEPTH) 1548 C BANDED EQUATION SOLVER -- CHOLES -- 1549 C IM=MAIN BAND WIDTH 1550 C IV=VERT BAND WIDTH 1551 C IA=AUG BAND WIDTH 1552 C C= CONSTANTS FOR THE NORMALS 1553 C T = NORMALS COEFFICIENTS 1554 C X = SOLUTION VECTOR 1555 C THIS SUBROUTINE REQUIRES FUNCTION JSUB 1556 C START OF THE FORWARD SOLUTION BY CHOLESKY METHOD 1557 IW=03 IW=06 0131 C**********CHECK DIAGONAL TERM FOR ZERO OR NEGATIVE 1559 C**********ALSO FIND MAX AND MIN VLUE OF DIAGONAL TERM 1560 IFIRST = 0 1561 700 CONTINUE 1562 XMIN = T(1,1) 1563 XMAX = XMIN 1564 DO 560 I=1,N 1565 II = JSUB(I,I,IM,IV,IA,N) 1566 IF(T(I,II)) 610,610,620 1567 610 WRITE(IW,515) I,I,T(I,II) 1568 515 FORMAT( 53H0ERROR 51 - NEGATIVE OR ZERO DIAGONAL TERM IN NORMALS/ 1569 1 3H T(,I3,1H,,I3,3H) =, E20.10) 1570 620 IF(XMIN-T(I,II)) 530,530,540 1571 540 MIN = I 1572 XMIN = T(I,II) 1573 GO TO 560 1574 530 IF(XMAX-T(I,II)) 550,550,560 1575 550 XMAX = T(I,II) 1576 MAX = I 1577 560 CONTINUE 1578 IF(IFIRST) 710,710,720 1579 710 WRITE (IW,715) 1580 715 FORMAT (26H0DIAGONAL TERMS OF NORMALS) 1581 GO TO 730 1582 720 WRITE (IW,725) 1583 725 FORMAT (34H0DIAGONAL TERMS OF NORMALS REDUCED) 1584 730 CONTINUE 1585 WRITE (IW,565) MIN,MIN,XMIN, MAX,MAX,XMAX 1586 565 FORMAT ( 1587 2 9H MIN = T(,I3,1H,,I3,3H) =, E20.10/ 1588 1 9H MAX = T(,I3,1H,,I3,3H) =, E20.10/) 1589 IF(XMIN) 570,570,580 1590 570 STOP 1591 580 CONTINUE 1592 IF(IFIRST) 740,740,800 1593 740 CONTINUE 1594 C**********END OF DIAGONAL TERM TEST 1595 T(1,1) = SQRT( T(1,1)) 1596 C(1) = C(1)/T(1,1) 1597 MU = IM+IV+IA 1598 DO 5 J=2,MU 1599 5 T(1,J) = T(1,J)/T(1,1) 1600 DO 100 I =2,N 1601 IF(I-(N-IV))20,20,10 1602 10 LLM = 1 1603 GO TO 30 1604 20 LLM = I-IM +1 1605 IF(LLM) 21,21,30 1606 21 LLM = 1 1607 30 SUMI = 0. 1608 ILLM = LLM 1609 SUMK=0. 1610 MU = I-1 1611 IF(MU-LLM)41,31,31 1612 31 DO 40 M = LLM,MU 1613 II = JSUB(M,I,IM,IV,IA,N) 1614 SUMI= SUMI+T(M,II)*T(M,II) 1615 40 SUMK = SUMK + T(M,II)*C(M) 1616 41 II = JSUB(I,I,IM,IV,IA,N) 1617 IF(T(I,II)-SUMI) 500,500,520 1618 500 WRITE(IW,510)I, I,T(I,II),SUMI 1619 510 FORMAT( 50H ERROR 50 - NEGATIVE ARGUMENT ON DIAGNAL SUMMATION / 1620 1 3H T(,I3,1H,, I3, 3H) =, E20.10/ 4H SUM, 8X, 1H=, E20.10) 1621 GO TO 570 520 CONTINUE 1622 T(I,II)=SQRT(T(I,II)-SUMI) 1623 C(I) = (C(I)-SUMK) /T(I,II) 1624 J=I+1 1625 IF(J -(N-IV))90,90,60 1626 90 IF((J-I)-IM)50,95,95 1627 50 LLM = J-IM +1 1628 IF(LLM)51,51,70 1629 51 LLM = 1 1630 GO TO 70 1631 60 LLM = ILLM 1632 70 SUMI = 0. 1633 IF(J-(N+IA))71,71,100 1634 71 IF(MU-LLM)81,72,72 1635 72 DO 80 M=LLM,MU 1636 IJ=JSUB(M,I,IM,IV,IA,N) 1637 JJ = JSUB(M,J,IM,IV,IA,N) 1638 80 SUMI = SUMI + T(M,IJ)*T(M,JJ) 1639 81 IJ = JSUB(I,J,IM,IV,IA,N) 1640 T(I,IJ) = (T(I,IJ)-SUMI)/T(I,II) 1641 J = J+1 1642 IF(J-(N-IV))90,90,60 1643 95 J = N-IV+ 1 1644 GO TO 60 1645 100 CONTINUE 1646 C BACK SOLUTION STARTS NOW 1647 NN = JSUB(N,N,IM,IV,IA,N) 1648 X(N) = C(N)/T(N,NN) 1649 I = N-1 1650 198 MLL = I+1 1651 M = MLL 1652 SUMI = 0. 1653 199 IF(M-(N-IV))201,201,202 1654 201 IF(M-I-(IM-1))202,202,200 1655 202 MI = JSUB(I,M,IM,IV,IA,N) 1656 SUMI = SUMI + T(I,MI)*X(M) 1657 GO TO 205 1658 200 M= N - IV 1659 205 M = M+1 1660 IF(M-N)199,199,203 1661 203 II = JSUB(I,I,IM,IV,IA,N) 1662 X(I) = (C(I)-SUMI)/T(I,II) 1663 I = I - 1 1664 IF (I) 204,204,198 1665 C NOW CALCULATE THE INVERSE ELEMENTS 1666 204 IF(IA) 400,400,300 1667 300 DO 320 I = 1, IA 1668 KL = IA + 1 1669 J=IM+IV+IA 1670 IN = I + IV+IM 1671 309 SUMI = 0.0 1672 KL = KL -1 1673 DO 310 K = KL , N 1674 310 SUMI = SUMI + T(K,IN)*T(K,J) 1675 T(I,J) = SUMI 1676 J = J-1 1677 IF(J-(IM+IV+I))320,309,309 1678 320 CONTINUE 1679 400 CONTINUE 1680 IFIRST = 2 1681 GO TO 700 1682 800 CONTINUE 1683 RETURN 1684 END 1685 FUNCTION JSUB(I,J,IM,IV,IA,N) 1688 IF(J-(N-IV)) 60,60,30 1689 30 JSUB = J - N + IV+ IM 1690 RETURN 1691 60 JSUB = (J-I) + 1 1692 RETURN 1693 END 1694 SUBROUTINE CURVE (IFIRST,T,NUMBER,NUMINV, RSL,RCL,X,NUMFR, 1697 1NUMTO,IA,IFIX,IG,ICODE,IV,IASTAB,IASTRO,IW,IPASS,IR,SD,IM, NSTA, 1698 2 NDEPTH, NWIDTH, IHEAD, LASTRO, LINV, RLONG, RLAT) 1699 DIMENSION T(NDEPTH,NWIDTH), NUMINV(LINV,2), NUMBER(NSTA,2), 1700 1 RSL(NSTA), RCL(NSTA), X(NDEPTH), NUMFR(2), NUMTO(2), IASTAB( 1701 2 LASTRO,2), IHEAD(20) 1702 DIMENSION RLONG(NSTA), RLAT(NSTA) 1703 DOUBLE PRECISION RLONG, RLAT, RSL, RCL 1704 IF(IFIRST) 20,30,20 1705 30 CONTINUE 1706 WRITE (IW,1115) IHEAD 1707 1115 FORMAT (1H1,15X,40HGEODETIC SURVEY OF CANADA - PROGRAM GALS// 1708 1 10X, A3, 19A4 ) 1709 WRITE (IW,05) SD 1710 05 FORMAT ( 54H0ELLIPSES BASED ON STANDARD DEVIATION OF UNIT WEIGHT = 1711 1, F7.3, 4H PPM//2X,2(7HSTATION,4X),2(2X,4HSEMI,7X),5HAZ OF, 1712 2 5X,6HERR AT, 4X, 5HRATIO, 5X, 2HAZ,6X, 4HDIST/ 1713 2 4X, 4HFROM, 7X, 2HTO, 6X, 10HMAJOR 1714 3AXIS,3X, 10HMINOR AXIS, 3X, 8HMAJ AXIS, 3X, 6H45 DEG, 3X, 1715 4 7HMIN/MAJ, 2(3X, 4HERR.,2X)) 1716 IFIRST = 1 1717 20 CONTINUE 1718 LLIM = IM+IV 1719 C ZERO X ARRAY 1720 DO 10 I=1,10 1721 10 X(I)=0 1722 C ZERO RADIUS OF CURVATURE 1723 XMP = 0 1724 PP = 0 1725 PQ = 0 1726 XMQ = 0 1727 IERROR = 0 1728 IENTER = 2 1729 CALL IGNOS (IENTER,IFR,ITO,NUMBER,NUMFR,NUMTO,IG,IERROR, 1730 1 IR,IW,IPASS,IERR,IASTAB,IASTRO,ICODE, NSTA,LASTRO,LINV) 1731 IF(IERROR) 100,15,100 1732 15 IF(ITO - IFR) 16,16,17 1733 16 INDTEM = ITO 1734 ITO = IFR 1735 IFR = INDTEM 1736 DO 18 I = 1,2 1737 INDTEM = NUMTO(I) 1738 NUMTO(I) = NUMFR(I) 1739 18 NUMFR(I) = INDTEM 1740 17 INDTEM = IA/2 1741 CALL SEARCH (IND,NUMINV,NUMFR,INDTEM, LINV ) 1742 CALL SEARCH (KEY,NUMINV,NUMTO,INDTEM,LINV ) 1743 IF(ITO-IFIX) 41,41,42 1744 41 IERROR = 60 1745 WRITE (IW,45) IERROR,(NUMFR(I),I=1,2),(NUMTO(I),I=1,2) 1746 45 FORMAT (6H ERROR,I3,17H - BOTH STATIONS ,2A4, 5H AND ,2A4, 1747 1 10H ARE FIXED) 1748 RETURN 1749 42 IF(IND) 50,60,50 1750 60 IF(IFR-IFIX) 61,61,62 1751 61 IF(KEY) 70,63,70 1752 63 IFR = ITO 1753 62 IERROR = 61 1754 WRITE(IW,65) IERROR,(NUMBER(IFR,I),I=1,2) 1755 65 FORMAT (6H ERROR,I3,35H - NO INVERSE ELEMENTS FOR STATION ,2A4) 1756 RETURN 1757 50 X(1) = T(IND*2-1,LLIM+IND*2-1) 1758 X(2) = T(IND*2-1,LLIM+IND*2) 1759 X(3) = T(IND*2,LLIM+IND*2) 1760 CALL RADIUS (XMP,PP,RSL(IFR),RCL(IFR)) 1761 GO TO 61 1762 70 X(8) = T(KEY*2-1,LLIM+KEY*2-1) 1763 X(9) = T(KEY*2-1,LLIM+KEY*2) 1764 X(10) = T(KEY*2,LLIM+KEY*2) 1765 CALL RADIUS (XMQ,PQ,RSL(ITO),RCL(ITO)) 1766 IF(IND) 80,90,80 1767 80 X(4) = T(IND*2-1,LLIM+KEY*2-1) 1768 X(5) = T(IND*2-1,LLIM+KEY*2) 1769 X(6) = T(IND*2,LLIM+KEY*2-1) 1770 X(7) = T(IND*2,LLIM+KEY*2) 1771 90 CALL ERCURV(XMP,XMQ,PP,PQ,X,NUMBER,IFR,ITO,IFIRST,IW,SD, 1772 1 NSTA, NDEPTH, RLONG, RLAT, RSL, RCL) 1773 100 RETURN 1774 END 1775 SUBROUTINE ERCURV ( XMP,XMQ,PP,PQ,X,NUMBER,IFR,ITO,IFIRST,IW,SD, 1778 1 NSTA, NDEPTH, RLONG, RLAT, RSL, RCL) 1779 DIMENSION X(NDEPTH), NUMBER(NSTA,2) 1780 DIMENSION RLONG(NSTA), RLAT(NSTA), RSL(NSTA), RCL(NSTA) 1781 DOUBLE PRECISION RLONG,RLAT,RSL,RCL, AZI, TWOTAN,SA12,CA12,SA21, 1782 1CA21,R1,W1,R2,W2,DIST 1783 PI = 3.141592654 1784 TWOPI = 3.141592654*2. 1785 C STANDARD DEVIATION IS IN PPM 1786 20 SO2 = SD*SD*1.E-12 1787 UA = SO2*(XMP*XMP*X(1) - 2*XMP*XMQ*X(4) + XMQ*XMQ*X(8)) 1788 UB = SO2*( PP*XMP*X(2) - XMP*PQ*X(5) - PP*XMQ*X(6) +XMQ*PQ*X(9)) 1789 UC = SO2*(PP*PP*X(3) -2*PP*PQ*X(7) + PQ*PQ*X(10)) 1790 BT = UA*UC - UB*UB 1791 AT = UA + UC 1792 SURD = SQRT(AT*AT - 4.*BT) 1793 AL = (AT + SURD)/2. 1794 BL = (AT - SURD)/2. 1795 S45 = SQRT(0.5*(AL+BL)) 1796 W = (ATAN((AL-UA)/UB)) 1797 CALL GEOINV(RLONG(IFR ),RLONG(ITO ),RSL(IFR ),RCL(IFR ),RSL(ITO 1798 1),RCL(ITO ),R1,W1,R2,W2,SA12,CA12,SA21,CA21,DIST) 1799 AZI =TWOTAN(SA12,CA12) 1800 ALFA = TWOPI - AZI - W 1801 SDIST = SQRT(AL*COS(ALFA)**2 + BL*SIN(ALFA)**2) 1802 SDIST = (SDIST/DIST) * 1.E06 1803 ALFA = ALFA - PI/2. 1804 SAZI = SQRT(AL*COS(ALFA)**2 + BL*SIN(ALFA)**2) 1805 SAZI = SAZI / (DIST * 4.85E-06) 1806 W = W*(-57.295779) 1807 SAL = SQRT(AL) 1808 SBL = SQRT(BL) 1809 RATIO = SBL/SAL 1810 WRITE (IW, 15) (NUMBER(IFR,I),I=1,2),(NUMBER(ITO,I),I=1,2), SAL, 1811 1 SBL,W,S45,RATIO, SAZI, SDIST 1812 15 FORMAT (1X,2(2A4,3X),F7.3,F13.3,F12.2,F10.2,F8.2, F10.2,F9.2) 1813 RETURN 1814 END 1815 SUBROUTINE RADIUS (RM,RP,SINPHI,COSPHI) 1818 DOUBLE PRECISION SINPHI,COSPHI,E2 1819 A = 6378206.4 1820 E2= .6768657997D-2 1821 D = 1.-E2*SINPHI*SINPHI 1822 RM = A*(1.-E2)/(D*SQRT(D)) 1823 RP = A*COSPHI/SQRT(D) 1824 RETURN 1825 END 1826