C *********************************************************** 00010025 C *********************************************************** 00020025 C ** ** 00030025 C ** THIS FORTRAN ROUTINE WAS MODIFIED BY J.R.ADAMS OF ** 00040025 C ** MCELHANNEY SURVEYING AND ENGINEERING LTD. TO RUN ** 00050025 C ** ON AN IBM COMPUTER USING FORTRAN G1. THE WORK WAS ** 00060025 C ** DONE FOR THE INDONESIAN GOVERNMENT DURING 1980 AND ** 00070025 C ** 1981. JOB # 083506 ** 00080025 C ** ** 00090025 C *********************************************************** 00100025 C *********************************************************** 00110025 C 00120025 C 00130025 C 00140021 C ******************************************************************00150021 C 00160021 C 00170021 C PROGRAM G * E * O * D * O * P 00180021 C 00190021 C *** V E R S I O N I I I (OCT. 74) *** 00200025 C BY J. KOUBA 00210011 C 00220012 C 00230018 C ******************************************************************00240019 C 00250020 C 00260021 C 00270021 C MXPS MAX NO OF ALLOWED POINTS 00280021 C MXDM MAX X DIMENSION I.E MXPS*3 00290021 C MYDM- Y DIMENSION MYDM= 3+ NOBI*MXPS 00300021 C NOBI -NO. OF BIASES IN X 00310021 C MXDOP MAX NO OF DOPLERS PER POINT 00320021 C MROW= MXDOP * MXPS 00330021 C 00340021 DOUBLE PRECISION A1,ESQ,PI,CMIN,CXYZ,ORBIT,DT,FS 00350021 DOUBLE PRECISION RHO,WE,RADG 00360021 C 00370021 DIMENSION R( 35707 ) , INT ( 618 ) 00380022 DIMENSION HEADER(9) 00390021 DIMENSION CHS(30),DXS(3) 00400021 DIMENSION DT(4),ORBIT(7),CXYZ(8,3) 00410021 C 00420021 C 00430021 C WHERE N1 IS GREATER OF MXDM OR MYDM 00440021 C 00450021 C INTDM= MXPS*MXDOP +6 MXPS+ MXDOP 00460021 C 00470021 C VARIABLE DIMENSION CHANGE THE FOLLOWING STATMENT IF REQUERED 00480021 C 00490021 C 00500021 C COMMON STATMENT BLOCKS 00510021 C BLANK COMM VARIABLE DIMENSION 00520021 C OPT COMM OPTIONS AND PROGRAM CONSTANTS 00530021 C COUNT COMM COUNT VARIABLES 00540021 C CHSQ COMM CHI-SQUARE ARRAY 00550021 C 00560021 C 00570021 REAL LO, HEADER 00580021 COMMON INT, R 00590021 COMMON /OPT/ ASG,HOR, SIGA, DXS,TE,XP99,RT,ALG,ACR,OUT,RHO, 00600021 1 FO, LO, WE, RADG, SGMT, ITTG, ITRIM, IWRT,IOUT,INPT, NROW 00610021 COMMON /COUNT/ NGRPTS, MACOL, MCCOL, MROW, IPASSC, IACCC, 00620021 1 IPASS, NREJT, NDFRJ , IELEV ,IVRJ , MDFM ,NPTS , ISEQ , IACC 00630021 COMMON /DBL/ A1,ESQ,PI ,CMIN,CXYZ 00640021 COMMON /CHI/ CHS 00650021 COMMON /SAT/ ORBIT,FS,DT,LOKR,IMIN,IST,NSTA,NDP,BIT,IT2,CPA 00660021 MAXDR = 35800 00670021 C 00680021 C 00690021 C PROGRAM CONSTANTS 00700021 C SGMT - A-PRIORI VAR. FACTOR 00710021 SGMT= 1. 00720021 C TROPCON - TRP. REFRACTION SCALING CONSTRAINT( 1 SIGMA) IN PER CENT00730021 TROPCO = 10. 00740021 C NOBI - NO. OF BIASES IN Y PARAMETERS 00750021 C MAX NO OF STA MXPS AND DOPPLERS CHANGE IF NECCESSARY 00760021 NOBI= 3 00770021 MXPS= 15 00780021 C SMXPS - MAXIMUM NO OF STATION OBSERVING SIMULTANEOUSLY A PASS 00790021 SMXPS= 15 00800021 MXDOP= 32 00810021 C 00820021 C CALCULATE THE REST 00830021 MXDM= MXPS* 3 00840021 MYDM= 3 + NOBI*SMXPS 00850021 MXROW= SMXPS* MXDOP 00860021 MXDOP2= MXDOP + 1 00870021 N1= MYDM 00880021 IF(N1.LT.MXDM) N1=MXDM 00890021 MYM=3+NOBI 00900021 N2= 7*MXDOP2+ MXDM +MXDM**2+MYDM**2 + MYDM 00910021 N4= (MXDOP**2 +MXDOP)/2 00920021 IF(N2.LT.N4) N2= N4 00930021 N3= MXROW + N1**2 + MYDM + N1*MYDM 00940021 N4= MXDOP**2 00950021 IF(N3.LT.N4) N3= N4 00960021 FO=1.92E 06 00970021 C 00980021 C INPT- INPUT (PREDOP OUTPUT FILE) 00990021 INPT= 2 01000021 C IOUT-OUTPUT UNIT NUMBER (PRINTER) 01010000 IOUT=6 01020001 1 IPASS=0 01030002 C O P T I O N C A R D READ IN 01040003 C 01050004 CALL OPTION (1,IEND,MRCN) 01060005 C 01070006 C TEST FOR EOF (IEND=1) 01080007 IF(IWRT.NE.10) REWIND 17 0109008 IF(IEND.EQ.1) GO TO 16 01100009 C 01110010 C APPRIORI SIGM SGMT CYCLES SQR,D 01120011 IPASSC=0 01130012 IACCC=0 01140013 IACC=0 01150014 MDFM=0 01160015 MCCOL= MYDM 01170016 MROW= MXROW 01180017 MACOL= MXDM 01190018 NPTS=0 01200019 NGRPTS=0 01210020 NJET=0 01220021 8 CONTINUE 01230021 C INITIALIZE ARRAYS 01240021 JH = MXDM**2 +23*MXPS+(4+MYM)*MXROW +2*MXPS*MXDOP+ MXDOP2 01250021 1+MXDM+ N2+ N3 01260021 DO 10 I=1 ,JH 01270021 10 R( I)=0. 01280021 WEIGHT= 625./ SGMT 01290021 PTRP= SGMT/TROPCO **2 01300021 J= MXDOP 01310021 K= 2*MXDOP 01320021 INTDM= MXPS*MXDOP +6*MXPS + MXDOP 01330021 DO 15 I=1,INTDM 01340021 15 INT(I)=0. 01350021 C 01360021 C CALCULATE ARRAY POSITIONS WITHIN R AND INT 01370021 C 'IDP' MULTP FACTOR FOR DOUBLE PRECISION 01380021 IDP = 2 01390021 IX = 1 +MXDM**2* IDP 01400021 IF = IX+ 6*MXPS * IDP 01410021 IG = IF + 5*MXPS 01420021 IS = IG + 3*MXPS * IDP 01430021 IA = IS + 3*MXPS 01440021 IC = IA + 3*MXROW * IDP 01450021 IW = IC + (3+NOBI)*MXROW * IDP 01460021 ID = IW + MXROW * IDP 01470021 IST= ID+ MXPS*MXDOP*IDP 01480021 IP =IST+ MXPS*MXDOP *IDP 01490021 ILK= IP + MXPS 01500021 IDM= ILK + MXPS 01510021 IEL=IDM+ 3*MXPS 01520021 IEX=IEL + MXDOP+1 01530021 IU= IEX +MXPS 01540021 ISK= IU + MXDM * IDP 01550021 IXS=ISK+(MXDOP+1) * IDP 01560021 IXX=IXS + 6*(MXDOP+1) * IDP 01570021 IPY=IXX+ MXDM 01580021 IYY= IPY + MYDM**2 * IDP 01590021 IAN=IYY+ MYDM * IDP 01600021 IV = IAN + MXDM**2 * IDP 01610021 IAM=IV + MXROW * IDP 01620021 IEND = IAM +((N1+1) * MYDM + N1*N1) * IDP 01630021 IF(IEND.GT.MAXDR) GO TO 16 01640021 C 01650021 KNG= 1+ MXPS*MXDOP 01660021 KNT=KNG+ MXDOP 01670021 KND=KNT+ MXPS 01680021 KNM=KND+ MXPS 01690021 KNR= KNM +MXPS *2 01700021 KNS= KNR + MXPS 01710021 KMG= KNS + MXPS 01720021 C 01730021 CALL CONTCD( MXPS, MXDM,R(IG), R(1),R(IX),VPV, R(IF), INT(KNT), 01740021 1 INT(KNM), IWRT, R(IS), R(IP), INT(KND), INT(KMG),SGMT) 01750021 C 01760021 IPASSC= IPASS 01770021 IACCC=IPASS- NDFRJ - NREJT 01780021 C 01790021 C 01800021 C READ IN FILE HEADER ON TAPE 2 AND CHEK FOR END OF FILE 01810021 C 01820021 READ(INPT,END=19) HEADER 01830021 17 WRITE(IOUT,18) HEADER 01840021 IF(IWRT.EQ.10) WRITE(17) HEADER 0185021 18 FORMAT(1H1,//,20H DATA SET FIGURE= ,2A4,4X,16HCHEB. FIT ORDER=, 01860021 1I3,5X,7HDATE = ,3A4,1X,16HNO. OF STATIONS=,I3,2X,15HEPHEM.SOURCE =01870021 1 ,2A4) 01880021 4 CONTINUE 01890021 C 01900021 CALL GEODOP ( MXDOP,MXDM,MXPS,MYDM, MXDOP2, MXROW, N1, R(1), R(IX)01910021 1, R(IF), R(IA), R(IC), R(IV),R(IAN), R(IW), R(IXX),R(IPY),R(IYY), 01920021 2 R(IAM),R(ISK), R(ID), R(IST), R(IEL), INT(1), INT(KNT), INT(KNM)01930021 3, R(IP), PTRP, R(IXS), INT(KND), R(IDM), R(IG), MYM, R(ILK), INT( 01940021 4KNR), R(IEX), R(IU), INT(KNS),INT(KNG),VPV, INT(KMG), MRCN) 01950021 C WRITE OPTIONS 01960021 C WRITE OPTIONS ON SUMMARY PAGE 01970021 CALL OPTION (0,IEND, MRCN) 01980021 C 01990000 WRITE(IOUT,18) HEADER 02000001 GO TO 20 02010002 C SUPPRESS STATION SUMMARY (IF APPLICABLE) 02020003 19 IF(ITRIM.LE.5) ITRIM= 6 02030004 20 CONTINUE 02040005 J2 = IAM + N1 * MYDM 02050006 JJ= JH - IF- MYDM 02060007 CALL SUMMAR ( MXDM, MXPS, R(1), R(IX), VPV, R(IF), R(IG), R(IU), 02070008 1R(IP), INT(KNT), R(IAN), INT(KNM), R(IEL), R(ISK), R(IPY), R(J2), 02080009 2 R(IS), INT(KNR), INT(KND), INT(KNS), R(ILK),JJ) 02090010 16 CONTINUE 02100012 STOP 02110013 END 02120014 BLOCK DATA 02130000 C 02140021 COMMON /OPT/ ASG,HOR, SIGA, DXS,TE,XP99,RT,ALG,ACR,OUT,RHO, 02150021 1 FO, LO, WE, RADG, SGMT, ITTG, ITRIM, IWRT,IOUT,INPT, NROW 02160021 COMMON /DBL/ A1,ESQ,PI,CMIN,CXYZ 02170021 COMMON /COUNT/ NGRPTS, MACOL, MCCOL, MROW, IPASSC, IACCC, 02180021 1 IPASS, NREJT, NDFRJ , IELEV ,IVRJ , MDFM ,NPTS , ISEQ , IACC 02190021 COMMON /CHI/ CHS 02200021 COMMON /SAT/ ORBIT,FS,DT,LOKR,IMIN,IST, NSTA,NDP,BIT,IT2,CPA 02210021 DOUBLE PRECISION A1, ESQ, PI , CMIN , CXYZ , ORBIT , DT , FS 02220005 DOUBLE PRECISION RHO , WE , RADG 02230021 C 02240006 DIMENSION CHS(30),DXS(3) 02250012 REAL L 02260020 DATA CHS/7.88,10.6,12.8,14.9,16.8,18.6,20.3,22.,23.6,25.2,26.8,28.02270000 13,29.8,31.3,32.8,34.3,35.7,37.1,38.6,40.,41.4,42.8,44.1,45.5,46.9,02280010 248.4,49.8,51.1,52.4,53.7/ 02290020 DATA RHO,PI,WE,RADG/2.06264806247D 05,3.141592653D0 , 4.3752691D-002300021 13 ,.1745329252D-01/ 02310021 DIMENSION DT(4),ORBIT(7),CXYZ(8,3) 02320021 DATA DT/.53678625D 00,.4601025D 00,.53678625D 00,.4663289583D00/ 02330021 DATA CMIN/1.798755D10 / 02340000 END 02350010 SUBROUTINE GEODOP (MXDOP, MXDM, MXPS ,MYDM, MXDOP2, MXROW, N1, PX,02360012 1 X, FDI, A, C, V, ANORM,W, XX, PY, YY, AMTR1, SK, DOP,ISHRT,02370013 2 ELEV, KDPSI, NTEST, NAME, PD, PTRP, XS, NDLY, DMET, GEOCO, NBI, 02380014 3LOK, NGR, ELMAX, U, NSEW, NG, VPV, MGR, MRCN) 02390015 C *********************************************************** 02400025 C *********************************************************** 02410025 C ** ** 02420025 C ** THIS FORTRAN ROUTINE WAS MODIFIED BY J.R.ADAMS OF ** 02430025 C ** MCELHANNEY SURVEYING AND ENGINEERING LTD. TO RUN ** 02440025 C ** ON AN IBM COMPUTER USING FORTRAN G1. THE WORK WAS ** 02450025 C ** DONE FOR THE INDONESIAN GOVERNMENT DURING 1980 AND ** 02460025 C ** 1981. JOB # 083506 ** 02470025 C ** ** 02480025 C *********************************************************** 02490025 C *********************************************************** 02500025 C 02510025 C 02520025 C MAIN PROGRAM SUBROUTINE 02530016 C 02540017 DOUBLE PRECISION A,AMTR1,ANORM,C,PX,PY,U,V,W,YY,DT,ISHRT 02550018 DOUBLE PRECISION A1,ESQ,GEOCO,PI,X,XS,DOP,SK,CMIN,ORBIT,FS 02560019 DOUBLE PRECISION RHO,WE,RADG,CXYZ 02570021 C 02580019 DIMENSION PX(MXDM,MXDM),X(6,MXPS),GEOCO(3,MXPS), FDI(5,MXPS), 02590010 1 A(MXROW,3 ), C(MXROW,NBI), W( 1), V( 1), ANORM(MXDM,MXDM), U(1 ),02600011 2 XX(1 ), PY(MYDM,MYDM), YY(1), SK(1), DOP(MXDOP,MXPS), ISHRT(MXDOP02610012 3,MXPS), NG(1), KDPSI(MXDOP,MXPS), NTEST( 1),NAME(2,1),XS(MXDOP2,6)02620013 C 02630014 DIMENSION PD(1), NDLY(1), NDESC( 4), DMET(3,MXPS), NGR(1), 02640015 1 CHS(30), LOK(1),DT(4),ORBIT( 7),ELMAX(1),DXS(3),CXYZ(8,3),NSEW(1)02650016 2,AMTR1(1), MGR(1) 02660017 DIMENSION ELEV(1) 02670017 C 02680018 C GEOCA LATITUDE, LONG. AND HEIGHT 02690019 C X X,Y,Z, AND TEPRORARY STORAGE AT ROW 4,5,6 02700020 C FDI FREQ. OFSET, CONSTRAINT,F-DRIFT, CONSTRAINT , PREVIOUS F.OFF02710021 C 02720021 C COMMON AREAS 02730021 C 02740021 REAL LOK 02750021 COMMON /OPT/ ASG,HOR, SIGA, DXS, TE,XP99,RT,ALG,ACR, OUT, 02760021 1RHO, FO,XLO, WE, RADG, SGMT, ITTG, ITRIM, IWRT,IOUT,INPT,NROW 02770021 COMMON /DBL/ A1,ESQ,PI,CMIN,CXYZ 02780021 C 02790021 COMMON /COUNT/ NGRPTS, MACOL, MCCOL, MROW, IPASSC, IACCC, 02800021 1 IPASS, NREJT, NDFRJ , IELEV ,IVRJ , MDFM ,NPTS , ISEQ , IACC 02810021 C 02820021 C 02830021 COMMON /SAT/ ORBIT,FS,DT,LOKR,IMIN,IST, NSTA,NDP,BIT,IT2,CPA 02840021 C 02850021 COMMON /CHI/ CHS 02860021 C 02870021 C INTERVAL LENGTH IN MIN AND C IN MIN 02880021 DABS(XX)=ABS(XX) 02890021 DCOS(XX)=COS(XX) 02900021 DSIN(XX)=SIN(XX) 02910021 DSQRT(XX)=SQRT(XX) 02920021 C 02930021 C 02940021 IACC= IACCC 02950021 LOKPR=0 02960021 IT2= MXDOP/8 02970021 BIT= 0.327708333E-03 02980021 DO 5 IA=1, MXPS 02990021 DO 1 I=1,MXDOP 03000021 L= I-((I-1)/IT2)*IT2 03010021 1 ISHRT(I,IA)= DT(L) 03020021 C MGR ACTUAL NO OF DOPPLERS READ IN 03030021 MGR(IA)= 0 03040021 NGR(IA)=0 03050021 ELMAX(IA)=0. 03060021 5 CONTINUE 03070021 XP100= 4.47 *DSQRT(SIGA) 03080021 7 CONTINUE 03090021 ISTOP=1 03100021 IZERO=1 03110021 DVPV=VPV 03120021 DVPV1=DVPV 03130021 C 03140021 CALL RINPT(MXPS, MXDOP,NGR,DOP,NDESC,ISHRT,IEND,NAME,KDPSI,NG, 03150021 1 DMET, LOK, NSEW, MGR, INJ , NTEST, LOKPR) 03160021 C 03170021 IF(IEND) 167,167,64 03180021 167 XLIM= XP100 03190021 168 CONTINUE 03200021 C ISTA- ACTUAL NO OF STATINS OBSERVED SIMULTANEOUSLY(.LE. MXPS) 03210021 ISTA=0 03220021 MRW=1 03230021 DMIN= (LOKR-LOKPR) /1440 03240021 IF(DMIN.GT.50)DMIN=0 03250021 DO 37 I=1, NGRPTS 03260021 DO 3 IB=1,3 03270021 3 X(IB+3,I)= X(IB,I) 03280021 C IGNORE STATION IF NTEST=10 03290021 C CHECK FOR MINIMUM DATA (2.5 MIN) 03300021 IMIND= (2.5/ISHRT(1,I) +.2) 03310021 IF(NGR(I).LT.IMIND) NGR(I)= 0 03320021 IF( NGR(I).EQ.0) GO TO 37 03330021 C UPDATE APPROXIMATE FREQ-CY FDI 03340021 FDI(1,I)= FDI(5,I) + FDI(3,I)*DMIN* .40E -01 03350021 C 03360021 CALL STCORD( ISHRT(1,I), MXDOP2,X(1,I), NDLY(I), FDI(1,I), SK, 03370021 1NG, NGR(I),KDPSI(1,I),ELMAX(I),ELEV,XS, LOK(I), GEOCO, DOP(1,I),V,03380021 1 NAME(1,I),MGR(I), NTEST(I) ) 03390021 C 03400021 IF( ELMAX(I).GT. TE) GO TO 6 03410021 NGR(I)= 0. 03420021 IELEV= IELEV+1 03430021 WRITE(IOUT,151) ELMAX(I), NAME(1,I),NAME(2,I), TE, IELEV 03440021 GO TO 37 03450021 6 NA=MRW+NGR(I)-1 03460021 DO 9 IA=MRW,NA 03470021 W(IA)=0. 03480021 V(IA)=0. 03490021 DO 8 IB=1,NBI 03500000 8 C(IA,IB)=0. 03510001 DO 9 IB=1,3 03520002 9 A(IA,IB)=0. 03530003 C 03540004 CALL FORM(NBI, MXDOP2, MXROW, A, C, W, X(1,I), GEOCO(3,I), 03550005 1FDI(1,I),SK,DOP(1,I), KDPSI(1,I), ISHRT(1,I), NTEST(I), PD(I), LOK03560006 1PR,NAME(1,I),ELEV,NGR(I), V, DMET(1,I), XS, LOK(I), NDLY(I),MRW, 03570007 , MGR(I), NG) 03580008 C 03590009 ISTA= ISTA +1 03600010 37 CONTINUE 03610011 NDD=MXDOP 03620012 C RESET VPV + VPV1 FOR ITERATION 03630013 VPV=DVPV 03640014 VPV1=DVPV1 03650015 NBO= NBI-3 03660016 MCCOL= 3+ ISTA*NBO 03670017 MROW= MRW - 1 03680018 DO 10 IA=1,MACOL 03690019 U(IA)=0. 03700020 XX(IA)=0. 03710021 DO 10 IB=1,MACOL 03720021 10 ANORM(IA,IB)=0. 03730021 DO 11 I=1,MCCOL 03740021 V(I)=0. 03750021 DO 12 J=1,MCCOL 03760021 12 PY(I,J)=0. 03770021 11 PY(I,I)= PTRP 03780021 C 10 P.C. TR. REFRACTION CONSTRAIN 03790021 C CONSTRAINT TIME BIAS UNIT 1=100MICROSEC 03800021 K=4 03810021 DO 28 I=1,NGRPTS 03820021 IF(NGR(I).EQ.0) GO TO 28 03830021 PY( K,K) = SGMT/FDI(2,I)**2/3600. 03840021 PY( K+1,K+1)= SGMT/FDI(4,I)** 2 03850021 K=K+NBO 03860021 28 CONTINUE 03870021 I= 1 03880021 PY(I ,I )= SGMT/ACR**2 03890021 PY(I+1,I+1)= SGMT/ALG**2 03900021 PY(I+2,I+2)= SGMT/OUT**2 03910021 NDGF= MROW - ISTA 03920021 IF(NDGF.GT.0) GO TO 147 03930021 NDFRJ=NDFRJ+1 03940021 WRITE(IOUT,148) NDFRJ 03950021 GO TO 7 03960021 147 CONTINUE 03970021 C 03980021 J1= 1+ N1* MYDM 03990021 J2= J1+ N1**2 04000021 C P H A S E A D J U S T M E N T C O M P U T A T I O N 04010021 C 04020021 CALL ADJST( MXDM, MYDM, MXROW, PX, A, C, W, U, ANORM, V, XX,04030021 1PY,YY,VPV,VPV1,PD, AMTR1(J2),AMTR1,AMTR1(J1), N1,NGR,NBI,NTEST) 04040021 C 04050021 C DATA SCREENING 04060021 2 IZERO= 1 04070021 IR=1 04080021 C ITRIM= 10 SIMULATION MODE 04090021 IF(ITRIM.EQ.10) GO TO 22 04100021 DO 20 K=1,NGRPTS 04110021 L= NGR(K) 04120021 IF(NGR(K).LE.0) GO TO 20 04130021 C TEST RESIDUALS FOR OUTLAYERS 04140021 C 04150021 CALL REJECT(V(IR),KDPSI(1,K),NG,DOP(1,K),MGR(K),IVRJ,XLIM,NGR(K) 04160021 1 ,ISTOP, IZERO,NAME(1,K)) 04170021 C 04180021 20 IR= IR +L 04190021 IF(ISTOP.EQ.0) GO TO 22 04200021 ISTOP=0 04210021 XLIM= XP99 04220021 IF(IZERO.EQ.0) GO TO 168 04230021 GO TO 2 04240021 22 IF(IZERO.EQ.0.AND.XLIM.NE.XP100) GO TO 167 04250021 C 04260021 C 04270021 C UPDATE APPROX. (PREVIOUS) CART. COORDINATES 04280021 C FREQ. OFFSET 04290021 C 04300021 IFF=0 04310021 C 04320021 C END ITERATION LOOP 04330021 C 04340021 JROW=1 04350021 PRINT 96 04360021 ISTA=1 04370021 DO 47 JK=1,NGRPTS 04380021 IF(NTEST(JK).EQ.10) GO TO 47 04390021 MADOP= MGR(JK) 04400021 JT= 4+ NBO*(ISTA-1) 04410021 IF(ITRIM.EQ.10) GO TO 48 04420021 DO 40 ICR=1,3 04430021 40 X(ICR,JK)= X(ICR+3,JK) + XX(ICR-1+JROW) 04440021 JROW= JROW +3 04450021 48 IF(NGR(JK).NE.0) GO TO 41 04460021 IFF=IFF+1 04470021 GO TO 47 04480021 41 FDI(1,JK)= FDI(1,JK) +YY(JT)/60.000 04490021 FOFF= FDI(1,JK) 04500000 IDLY= NDLY(JK) +YY(JT +1)*100 04510001 J= JT +2 04520002 LLOK= LOK(JK) 04530003 CALL DIME(LLOK,ID,IH,IMN) 04540004 46 WRITE(IOUT,95)NAME(1,JK),IH,IMN,(DMET(JL,JK),JL=1,3),FOFF,IDLY,YY(04550005 1 J),ELMAX(JK), (X(JL,JK),JL=1,3), (KDPSI(JL,JK),JL=1,MADOP) 04560006 ISTA= ISTA +1 04570007 47 CONTINUE 04580008 C 04590009 IF(ITTG.EQ.5) GO TO 175 04600010 WRITE(IOUT,126) 04610011 DO 145 J=1,3 04620012 145 WRITE(IOUT,98)(XX(I*3-3+J),I=1,NGRPTS) 04630013 175 CONTINUE 04640014 C PRINT VECTOR V 04650015 C 04660016 DTVPV=VPV-DVPV 04670017 C COMPUTE NO. OF DEGREES OF FREEDOM 04680018 C 04690019 C STATISTICAL TESTING OF VPV (CHI-SQUARE TEST) 04700020 C 04710021 IDFM=MDFM 04720021 MDFM=MDFM+NDGF 04730021 IDFM=MDFM-IDFM 04740021 IF (IDFM.LE.0) GO TO 52 04750021 IF(ITRIM.EQ.10) GO TO 54 04760021 CHISQ= SIGA*(IDFM + SQRT(2.* IDFM)*XP99) 04770021 IF(IDFM.LE.30) CHISQ= SIGA* CHS(IDFM) 04780021 IF (DTVPV.LE.CHISQ) GO TO 54 04790021 52 MMMM=MDFM-IDFM 04800021 DO 53 JK=1,NGRPTS 04810021 DO 51 K=1,3 04820021 51 X(K,JK)= X(K+3,JK) 04830021 53 FDI(1,JK)= FDI(5,JK) 04840021 MDFM=MMMM 04850021 VPV=DVPV 04860021 VPV1=DVPV1 04870021 NREJT=NREJT+1 04880021 WRITE(IOUT,89) (YY(I),I=1,3),VPV,DTVPV,MDFM, SIGM 04890021 PRINT 90, NREJT 04900021 GO TO 7 04910021 54 CONTINUE 04920021 SIGM=1. 04930021 IF(ITRIM.EQ.10) GO TO 66 04940021 IF(MDFM) 66,66,67 04950021 67 SIGM= VPV/MDFM 04960021 66 CONTINUE 04970021 C 04980021 C ADJUSTMENT 04990021 IO= -N1 +J1 -1 05000021 DO 55 I=1, MACOL 05010021 IO= IO + N1 05020021 DO 55 J= 1, MACOL 05030021 JI= IO + J 05040021 55 PX(J,I)= AMTR1(JI) 05050021 IF(ITTG.EQ.5) GO TO 63 05060021 PRINT 117 05070021 DO 57 II= 1, MACOL 05080021 57 WRITE(IOUT,125) (PX(II,K),K=II,MACOL ) 05090021 C 05100021 C ESTIMATION OF VARIANCE COVARIANCE MATRIX 05110021 C 05120021 C WRITE OUT RESIDUALS 05130021 WRITE(IOUT,128) 05140021 NGRS=1 05150021 DO 50 I=1,NGRPTS 05160021 IF(NGR(I).EQ.0) GO TO 50 05170021 NGRE= NGRS+NGR(I)-1 05180021 WRITE(IOUT,92) (NAME(J,I),J=1,2) , (V(J),J=NGRS,NGRE) 05190021 NGRS=NGRE+ 1 05200021 50 CONTINUE 05210021 WRITE(IOUT,91) 05220021 63 CONTINUE 05230021 DO 60 JK=1,MACOL 05240021 DO 59 IK=1,MACOL 05250021 59 ANORM(JK,IK)=SIGM*ANORM(JK,IK) 05260021 IF(ITTG.NE.5) WRITE(IOUT,125) (ANORM(JK,IK),IK=JK,MACOL) 05270021 60 CONTINUE 05280021 62 SIGM=DSQRT(SIGM) 05290021 WRITE(IOUT,89) (YY(I),I=1,3),VPV,DTVPV,MDFM, SIGM 05300021 IACC= IACC +1 05310021 WRITE(IOUT,177) IACC 05320021 IF(IWRT.EQ.10) GO TO 69 05330021 IF(ITRIM.EQ.10) GO TO 7 05340021 KT= MACOL 05350021 WRITE(17)IPASS,IST, (LOK(K), (X(KX,K), KX=1,3), FDI(1,K), NGR(K),0536021 1 NSEW(K),ELMAX(K),K=1,NGRPTS), ((ANORM(I,J),J=I,KT),I=1,KT),SIGM ,05370021 2 MDFM 05380021 GO TO 7 05390021 C WRITE OUT CLEANED INPUT(TAPE5) ONTO TAPE7 05400021 69 WRITE(17)NDESC, IST, INJ, LOKR, IMIN, FS, NSTA, ORBIT, NROW, CXYZ,0541021 1 (MGR(I), I=1,NSTA) 05420021 DO 68 I=1,NSTA 05430021 NDD= MGR(I) 05440021 IT= NDD 05450021 IF(MRCN.EQ.1) IT=0 05460021 IF( NDD.EQ.0) GO TO 68 05470021 C WRITE(17)NAME(1,I),LOK(I), (DMET(K,I),K=1,3), NSEW(I), FDI(1,I), 0548021 WRITE(17)NAME(1,I),NAME(2,I),LOK(I),(DMET(K,I),K=1,3),NSEW(I), 0549021 1FDI(1,I),IT,( DOP(J,I),J=1,NDD), (ISHRT(J,I),J=1,IT) 05500021 68 CONTINUE 05510000 177 FORMAT('0', 5X,'--------------------------------------------------05520001 1---------- END OF ACCEPTED PASS NO=',I4,' ------------') 05530002 GO TO 7 05540003 64 ENDFILE 7 05550004 RETURN 05560005 C 05570006 151 FORMAT('0 ++++ PASS ELEVATION ',F4.1,' AT STATION ',2A4, ' LESS T05580007 1HAN ',F4.1,' DEG CRITER. -REJECTION NO ',I4) 05590008 148 FORMAT('0',' ++++ PASS REJECTED ON ZERO OR NEGATIVE D.F -REJECTION05600009 1 NO ',I4) 05610010 90 FORMAT('0',//,' ++++ PASS REJECTED AT 99 P.C. PROB.LEVEL -REJECTIO05620011 1N NO =',I4) 05630012 92 FORMAT(' ',2X,2A4,15F5.1/' ',10X,15F5.1) 05640026 91 FORMAT('0',//' ESTIMATED VARIANCE COVARIANCE MATRIX'/) 05650026 95 FORMAT(1H ,2X,A4,I3,I2,F8.0,2F5.1,F9.2,I6,2F6.1,2X,3F11.1,1X,8(1X,05660015 14I1)) 05670016 89 FORMAT(1H0,21H ORB BIASES(M) ACCR=,F6.1,5H ALG=,F6.1,5H OUT=,F6.105680017 1,9X,5H VPV=,F9.2,6H DVPV=,F7.2,5H DF =,I5,5X,4H SO=,F6.2) 05690018 96 FORMAT(1H0,56H ST.NO H MN PRESS TEMP PWP FREQ. DELAY TROP: 05700019 1ELEV,9X,1HX,10X,1HY,10X,1HZ,9X,14H USED DOPPLERS ) 05710020 98 FORMAT(15F9.2) 05720021 117 FORMAT('0',' PHASE WEIGHT COEFICIENT MATRIX PX ',/ ) 05730021 125 FORMAT(' ',3(3X,3E13.5)) 05740021 126 FORMAT('0',' S O L U T I O N V E C T O R X ') 05750021 128 FORMAT('0', / ,' VECTOR V') 05760021 END 05770021 SUBROUTINE CONTCD( N1, N2, GEO, PHWTM, X, VPV, FOFFS, NTEST, 05780024 1SNAME, IN, SHIFT, RCRMX, NDLY, MGR ,SGMT ) 05790024 C 05800024 C *********************************************************** 05810025 C *********************************************************** 05820025 C ** ** 05830025 C ** THIS FORTRAN ROUTINE WAS MODIFIED BY J.R.ADAMS OF ** 05840025 C ** MCELHANNEY SURVEYING AND ENGINEERING LTD. TO RUN ** 05850025 C ** ON AN IBM COMPUTER USING FORTRAN G1. THE WORK WAS ** 05860025 C ** DONE FOR THE INDONESIAN GOVERNMENT DURING 1980 AND ** 05870025 C ** 1981. JOB # 083506 ** 05880025 C ** ** 05890025 C *********************************************************** 05900025 C *********************************************************** 05910025 C 05920025 C 05930025 DOUBLE PRECISION DD,EE,GEO,PHWTM,TT,WEIGHT,X,CMIN,CXYZ 05940025 DOUBLE PRECISION SS,VV,WW 05950024 C 05960024 C WRITTEN BY P. LAWNIKANIS NOV 22/73 05970024 C REREREREMODIFIED FEB 12/74 05980024 C 05990024 C 06000024 C THIS ROUTINE WILL READ A .GEODOP. 06010024 C CONTINUATION FILE FROM CARD INPUT OR .TAPE8. - 06020024 C 06030024 C USE THE SUMMARY PAGES FROM THE LAST 06040024 C GOOD .GEODOP. RUN TO PREPARE THE DECK - 06050024 C 06060024 C 06070024 C A FORMAT GUIDE FOLLOWS ..... 06080024 C .I. = AN INTEGER DIGIT 06090024 C 06100024 C........1.........2.........3.........4.........5.........6.........7..06110024 C*1* I III III III III III IIIISF06120024 C NO.STNS/NO.PASSES/ RJT-99PC/ RJT-ZDF/ RJT-LOEL/ DOP-99PC/TOT NO.DF/ 06130024 85 FORMAT(1H0,2A4,2X,I2,F5.1,I5,F5.1,F8.1,2F5.1,2X,2(2I3,F7.3,1X), 06140024 1 F7.1,1X,3F5.1,1X,3F11.1) 06150024 C 06160024 C 06170024 C 06180024 COMMON /COUNT/ NSTNS,K,II,IJ,IK,IL,NPASS,N99PR,NZDFR, 06190024 $ NLELR,ND99R,NFREE,NNTNS,NNNNS,IACCP 06200024 COMMON /OPT/ AA,BB,CC, HH,HI,HJ,HK,HL,HN,OO,PP,QQ,SS, 06210024 1 ZZ, UU, VV, WW, XX, JJ, KF, LL, MM, NN, NNNN 06220024 COMMON /DBL/ DD,EE,TT , CMIN , CXYZ 06230024 C 06240024 DIMENSION FOFFS(5, N1), NTEST(1), PHWTM(N2,N2), RCRMX( 1), GEO(3,N06250024 11), SHIFT(3,N1), SNAME(2,1),X(6,N1),NDLY(1) 06260024 DIMENSION MGR(1) 06270024 DATA INPC /'INP?'/ , KAP8/'TAP8'/ , KAP9/'TAP9'/ 06280024 C 06290024 C INPUT CARD CONTINUATION FILE ..... 06300024 C ... NO. OF - STATIONS,PASSES,99 P.C. REJECTIONS,ZERO D.F.,LO EL. REJ06310024 C DOPPLERS REJECTED,TOTAL D.F.,AND VPV ... 06320024 C 06330024 C 06340024 C THIS SETTING OF 'PX' WAS REMOVED FROM MAINPGM 06350024 C 06360024 CCARD READER SELECT CODE 06370024 IREAD=5 06380024 WEIGHT = 1.0D4 / SGMT 06390024 DO 19 I = 1,N2 06400024 19 PHWTM(I,I) = WEIGHT 06410024 DO 28 I=1,N1 06420024 28 MGR(I)= N1 06430024 C 06440024 20 READ(IREAD,5) NSTNS,NPASS,N99PR,NZDFR,NLELR,ND99R,NFREE,VPV 06450024 5 FORMAT(8X,I2,5(7X,I3),4X,I6,F10.3) 06460024 N = INPC 06470024 C 06480024 C ... NAME,TEST,AND X Y Z FOR EACH STATION ..... 06490024 C 06500024 C NNTNS - NO. OF ACTIVE STATIONS (NTEST.NE.10) 06510024 NNTNS =0 06520024 DO 100 I=1,NSTNS 06530024 READ(IREAD,15) (SNAME(J,I),J=1,2), 06540024 + NTEST(I), RCRMX(I),NDLY(I), (FOFFS(J,I),J=1,4),06550024 1(SHIFT(J,I),J=1,3) 06560024 15 FORMAT(1X,2A4,1X,I5,F5.1,I5,F5.1,3F5.0,3F5.1) 06570024 READ(IREAD,16) (X(J,I),J=1,3) 06580024 16 FORMAT(3F10.1) 06590024 C 06600024 C CALCULATE LATITUDE,LONGDITUDE,AND ELEVATION FROM X Y Z ..... 06610024 C 06620024 IF(IN.EQ.8.OR.IN.EQ.9) GO TO 250 06630024 CALL CARGEO (X(1,I), X(2,I), X(3,I), GEO(1,I), GEO(2,I), GEO(3,I)06640024 1 ,YY, DD, EE) 06650024 250 CONTINUE 06660024 IF(NTEST(I).EQ.10) GO TO 27 06670024 C SET ACTIVE STATION INDEX MGR 06680024 NNTNS= NNTNS +1 06690024 MGR(NNTNS)= I 06700024 27 CONTINUE 06710024 C SET DEFAULT VALUES WEIGHTS OF DOPPLERS, DF, TIM. BIAS 06720024 IF(RCRMX(I).EQ.0.) RCRMX(I)= 1. 06730024 IF( FOFFS(2,I).EQ.0.) FOFFS(2,I)= 50. 06740024 IF( FOFFS(4,I).EQ.0.) FOFFS(4,I)= 1. 06750024 FOFFS(5,I)= FOFFS(1,I) 06760024 100 CONTINUE 06770024 K= NNTNS * 3 06780024 C 06790024 C ... UPPER TRIANGULAR ELEMENTS OF PHASE WEIGHT MATRIX ... 06800024 C 06810024 READ( IREAD,35,END=220) (PHWTM( I , I ) , I=1,K) 06820024 35 FORMAT(6(2X,G11.5)) 06830024 210 READ( IREAD,36,END=220) IRW , ICOL , NOEL 06840024 36 FORMAT(3I10) 06850024 240 IF((IRW*ICOL*NOEL).EQ.0) GO TO 220 06860024 IF(IRW.GE.K.OR.ICOL.GT.K) GO TO 220 06870024 C IRW - STARTING ROW 06880024 C ICOL - STARTING COL 06890024 C NOEL - NO. OF ELEMENTS TO BE READ IN BY ROWS (UPPER TRIANGLE) 06900024 230 ICNT=0 06910024 ICL= ICOL 06920024 KM1= K-1 06930024 DO 200 I=IRW,KM1 06940024 IF(ICNT.GE.NOEL) GO TO 210 06950024 KCL= ICL +NOEL -ICNT -1 06960024 IF(KCL.GT.K) KCL= K 06970024 READ( IREAD,35,END=220) (PHWTM(I,J),J=ICL,KCL) 06980024 205 CONTINUE 06990024 ICNT= ICNT + KCL- ICL +1 07000024 ICL= I +2 07010024 C 07020024 C ... SET LOWER TRIANGULAR ELEMENTS ... 07030024 C 07040024 200 CONTINUE 07050024 220 IF(PHWTM(1,1).EQ.0) PHWTM(1,1)=10000./XX 07060024 C 07070024 C CHECK INPUT MODE PARAMETR 07080024 PRINT 115 07090024 C 07100024 C 07110024 C IN - TAPE NO. = 9 (REWOUND) OR = 8(NOT REWOUND) 07120024 C 07130024 IF(IN -8) 40,10, 11 07140024 11 IF(IN.GT.9) GO TO 40 07150024 REWIND IN 07160024 N=KAP9 07170024 GO TO 12 07180024 10 N=KAP8 07190024 12 IOUT= IN 07200024 READ(IOUT,END=29) NNN,MACOL,(J,I=1,NNN), 07210024 + (SNAME(1,MGR(J)),SNAME(2,MGR(J)),(GEO(I,MGR(J)), 07220024 1 X(I,MGR(J) ), I=1,3), J=1,NNN), 07230024 2 ((PHWTM(J,I),I=J,MACOL),J=1,MACOL), NFREE, VPV, NPASS, NLELR, 07240024 2 NZDFR, N99PR, ND99R 07250024 GO TO 12 07260024 29 IF(IN.EQ.9) REWIND IN 07270024 C 07280024 C ... OUTPUT DATA TO PRINTER FOR VERIFICATION ... 07290024 C 07300024 40 PRINT 75,N,NPASS,NSTNS,N99PR,NZDFR,HK,NLELR,ND99R 07310024 PRINT 90 07320024 90 FORMAT(1H0,//,3X,45HST.NAME CODE WGHT DELY T.B. FREQ. STD. DRIFT07330024 1,3X, 8HLATITUDE,8X, 9HLONGITUDE,2X,29HHEIGHT LAT/LONG/HGT ANT SHIF07340024 1T,3X, 1HX,10X,1HY,10X,1HZ) 07350024 DO 300 I=1,NSTNS 07360024 C 07370024 C CALCULATE DEG. MIN. SEC. FOR PRINTOUT ..... 07380024 C 07390024 CALL QUICK( A, B, C, GEO(1,I)) 07400024 CALL QUICK( D, E, F, GEO(2,I)) 07410024 C 07420024 C MODIFY FREQ. OFFSET(S) IF READ IN ..... 07430024 C 07440024 SGN = 1.0 07450024 IF ( GEO(1,I) .NE. 0.D0 )SGN = GEO(1,I) / DABS(GEO(1,I)) 07460024 KK= A+.5 * SGN 07470024 L= B*SGN + .5 07480024 C = C * SGN 07490024 IF ( GEO(2,I) .NE. 0.D0 ) SGN = GEO(2,I) / DABS ( GEO(2,I)) 07500024 M= D+.5 * SGN 07510024 N= E * SGN +.5 07520024 F = F * SGN 07530024 H= GEO(3,I) 07540024 PRINT 85, (SNAME(J,I),J=1,2), 07550024 + NTEST(I), RCRMX(I), NDLY(I), FOFFS(4,I),(FOFFS07560024 1(J,I),J=1,3),KK,L,C,M,N,F,H,( SHIFT(J,I),J=1,3),(X(J,I),J=1,3) 07570024 300 CONTINUE 07580024 PRINT 95,VPV,NFREE 07590024 DO 206 I=1,K 07600024 DO 206 J=I,K 07610024 206 PHWTM(J,I)= PHWTM(I,J) 07620024 DO 400 I=1,K 07630024 PRINT 105,( PHWTM(I,J),J=I,K) 07640024 400 CONTINUE 07650024 IF( IN.EQ.1) GO TO 402 07660024 CALL TRINN(PHWTM, K, 2, N2) 07670024 DO 401 I=1,K 07680024 DO 401 J=I,K 07690024 401 PHWTM(J,I)= PHWTM(I,J) 07700024 402 CONTINUE 07710024 PRINT 25 07720024 PRINT 26 07730024 25 FORMAT(//,37H REMARKS- CODE=0 REC TIME FRAME(CBR)/15X,22H=1 SAT.T07740024 1IME FRAME( BR)/15X,22H=10 IGNORE THE STATION/11X,38HWGHT= WEIGHT( 07750024 11/(ST.DEV.IN COUNTS)**2)/11X, 07760024 129HDELY= RECEVER DELAY(1E-6 SEC)/11X,42HT.B.= TIMING BIAS CONST(1 07770024 2SIGMA) (1E-4SEC)) 07780024 26 FORMAT(1H , 10X,35HFREQ= FREQ. OFFSET FROM 400MHZ (HZ)/ 07790024 1 11X,36HSTD.= FR07800024 1EQ. OFFSET CON.(1SIGMA) (HZ)/11X,30HDRIFT= FREQ. DRIFT( 1E-10/DAY)07810024 3/2X,33HHEIGT, ANT. SHIFT, X,Y,Z (METERS)) 07820024 C RETURN TO .GEODOP. ..... 07830024 C 07840024 CALL CENTER(X, GEO, DD, EE, NSTNS, SHIFT,1.) 07850024 50 RETURN 07860024 C 07870024 75 FORMAT ( /,51X,'CONTINUATION DATA FROM ',A7,/,1H+,49X, 07880024 1////,' TOTAL NUMBER OF PROCESSED PASSES =',I4,/, 07890024 $ ' NUMBER OF UNKNOWN STATIONS =',I4,////07900024 $ /,' REJECTED ON 99 PER CENT =',I4,/, 07910024 $ ' REJECTED ON ZERO DEG FRDM =',I4,/, 07920024 + ' REJECTED ON ',F4.1,' DEG MAX EL =',I4,//, 07930024 $ ' DOPPLERS REJECTED ON 99 PER CENT =',I4, 07940024 1////) 07950024 95 FORMAT(//,' VPV = ',F10.3,/' DEGREES OF FREEDOM = ',I5,//, 07960024 1 ' WEIGHT COEFFICIENT MATRIX -' ) 07970024 105 FORMAT(1H ,(3(3X,3E12.5)) ) 07980024 115 FORMAT (1H1) 07990024 C 08000024 END 08010024 SUBROUTINE QUICK (A,B,C,D) 08020021 DOUBLE PRECISION D,E 08030021 C *********************************************************** 08040025 C *********************************************************** 08050025 C ** ** 08060025 C ** THIS FORTRAN ROUTINE WAS MODIFIED BY J.R.ADAMS OF ** 08070025 C ** MCELHANNEY SURVEYING AND ENGINEERING LTD. TO RUN ** 08080025 C ** ON AN IBM COMPUTER USING FORTRAN G1. THE WORK WAS ** 08090025 C ** DONE FOR THE INDONESIAN GOVERNMENT DURING 1980 AND ** 08100025 C ** 1981. JOB # 083506 ** 08110025 C ** ** 08120025 C *********************************************************** 08130025 C *********************************************************** 08140025 C 08150025 C 08160025 C 08170021 C THIS SUBROUTINE COVERTS A COORDINATE IN 08180021 C RADIANS TO ONE IN DEG. MIN. SEC. ..... 08190000 C 08200001 E = D * 57.2957795129D0 08210002 A = FLOAT(IDINT(E)) 08220003 E = ( E - A ) * 60. 08230004 B = FLOAT(IDINT(E)) 08240005 C = ( E - B ) * 60. 08250006 C 08260007 RETURN 08270008 END 08280009 SUBROUTINE CENTER(X, GEOCO, A1, ESQ, NST , SHFT, SGN) 08290015 C *********************************************************** 08300025 C *********************************************************** 08310025 C ** ** 08320025 C ** THIS FORTRAN ROUTINE WAS MODIFIED BY J.R.ADAMS OF ** 08330025 C ** MCELHANNEY SURVEYING AND ENGINEERING LTD. TO RUN ** 08340025 C ** ON AN IBM COMPUTER USING FORTRAN G1. THE WORK WAS ** 08350025 C ** DONE FOR THE INDONESIAN GOVERNMENT DURING 1980 AND ** 08360025 C ** 1981. JOB # 083506 ** 08370025 C ** ** 08380025 C *********************************************************** 08390025 C *********************************************************** 08400025 C 08410025 C SUBROUTINE CENTER 08420025 C USING 1ST ORDER DIFFERENTIAL FORMUALE (EXACT FOR HEIGHT) 08430025 C COMPUTES LAT, LOG AND HEIGHT FOR THE MONUMENT FROM ECC POSITION 08440025 C AND VICE VERSA 08450025 C SGN =1 FROM MONUMENT TO ANTENA =-1 THE OPPOSITE WAY 08460025 C 08470025 DOUBLE PRECISION A1,DCSF,DCSL,DSSL,DSSF,ESQ,GEOCO,RM,RN,W,X 08480021 DIMENSION X(6,NST), GEOCO(3, NST), SHFT(3,NST) 08490016 DO 1 I=1,NST 08500017 DCSF=DCOS(GEOCO(1,I)) 08510018 DCSL=DCOS(GEOCO(2,I)) 08520019 DSSL=DSIN(GEOCO(2,I)) 08530020 DSSF=DSIN(GEOCO(1,I)) 08540021 W=DSQRT(1.D0 - ESQ*DSSF*DSSF) 08550021 RN= A1/W + GEOCO(3,I) 08560021 RM= RN*( 1.0D0 - ESQ)/W**3 08570021 GEOCO(3,I)= GEOCO(3,I) + SGN* SHFT(3,I) 08580021 GEOCO(1,I)= GEOCO(1,I) + SGN* SHFT(1,I) /RM 08590021 GEOCO(2,I)= GEOCO(2,I) + SGN* SHFT(2,I) /(RN* DCSF) 08600021 C X - CENTERING 08610021 X(1,I)= X(1,I) +(-SHFT(1,I)*DSSF*DCSL-SHFT(2,I)*DSSL+SHFT(3,I)* 08620021 1DCSF*DCSL)* SGN 08630021 X(2,I)= X(2,I) +(-SHFT(1,I)* DSSF*DSSL + SHFT(2,I)*DCSL +SHFT(3,I)08640021 1 *DCSF*DSSL)*SGN 08650021 X(3,I)= X(3,I) +( SHFT(1,I)*DCSF +SHFT(3,I)*DSSF)* SGN 08660021 1 CONTINUE 08670021 RETURN 08680021 END 08690021 SUBROUTINE PPLLOT( N , COORD, NAME, NN, XYZF) 08700021 C *********************************************************** 08710025 C *********************************************************** 08720025 C ** ** 08730025 C ** THIS FORTRAN ROUTINE WAS MODIFIED BY J.R.ADAMS OF ** 08740025 C ** MCELHANNEY SURVEYING AND ENGINEERING LTD. TO RUN ** 08750025 C ** ON AN IBM COMPUTER USING FORTRAN G1. THE WORK WAS ** 08760025 C ** DONE FOR THE INDONESIAN GOVERNMENT DURING 1980 AND ** 08770025 C ** 1981. JOB # 083506 ** 08780025 C ** ** 08790025 C *********************************************************** 08800025 C *********************************************************** 08810025 C 08820025 C 08830025 DIMENSION COORD(5,NN), XYZF(5) 08840021 DATA RESOLC, RESOLF/ 12.,6./ 08850021 DIST(A,B,C)= SQRT(A**2+B**2+C**2) 08860021 C 08870021 C NN= MAXIMUM NO OF PASSES TO BE PLOTED 08880021 C E) G) FOR MXDOP= 8 08890021 C MXPS NN 08900021 C 1 50 08910021 C 2 120 08920021 C 3 250 08930021 C 5 600 08940021 C 08950021 DO 200 J=1,N 08960021 C PLOT X, Y, Z, DF 08970021 R= RESOLC 08980021 IF(J.EQ.1) R= RESOLF 08990021 ZT= FLOAT(INT(XYZF(J))) 09000021 XT= ZT-R 09010021 YT= ZT+R 09020021 CALL PRINTP(1,J,XT,YT,ZT, NAME,1,1) 09030021 DO 300 K=1,NN 09040021 300 CALL PRINTP(2,J,XT,XT,COORD(J,K),J,J,J) 09050021 CALL PRINTP(3,J,XT,XT,XT,J,J,J) 09060021 200 CONTINUE 09070021 RETURN 09080021 END 09090021 SUBROUTINE PRINTP (MODE,LABEL,XLOWR,XHIGH,VALUE,NAME1,NAME2,LINE)09100021 C *********************************************************** 09110025 C *********************************************************** 09120025 C ** ** 09130025 C ** THIS FORTRAN ROUTINE WAS MODIFIED BY J.R.ADAMS OF ** 09140025 C ** MCELHANNEY SURVEYING AND ENGINEERING LTD. TO RUN ** 09150025 C ** ON AN IBM COMPUTER USING FORTRAN G1. THE WORK WAS ** 09160025 C ** DONE FOR THE INDONESIAN GOVERNMENT DURING 1980 AND ** 09170025 C ** 1981. JOB # 083506 ** 09180025 C ** ** 09190025 C *********************************************************** 09200025 C *********************************************************** 09210025 C 09220025 C 09230025 DIMENSION CHAR(121),HEAD(2,5),IAXIS(13),NAME1(2) 09240021 DATA BLANK/4H / 09250021 DATA CHAR/121*1H / 09260021 DATA HEAD/'FREQ','.OFF' , '** D','X **' , '** D','Y **' , 09270021 ^ '** D','Z **' , '** D','R **' / 09280021 DATA IPLOT,JT/0,1/ 09290021 DATA BLA/' '/ , AST/'*'/ , HEX/'X'/ 09300021 C 09310021 C .MODE. = 1 ... START NEW PLOT WITH .LABEL.TH HEADING, 09320021 C WITH LIMITS .XLOWR. LE .VALUE. LE .XHIGH., 09330021 C STATION NAMES .NAME1. AND .NAME2., 09340021 C AND ABCISSA STARTING AT .LINE. . 09350021 C = 2 ... PLOT .VALUE. IF IN RANGE (*), OR (X) 09360021 C IF OUT-OF-RANGE, AND PRINT POINT NO. 09370021 C IF A MULTIPLE OF 5 . 09380021 C = 3 ... PRINT ORDINATE AXIS TO FINISH PLOT . 09390021 C 09400021 GO TO (10,20,30),MODE 09410021 C 09420021 10 IPLOT = IPLOT + 1 09430021 PRINT 5,IPLOT,(HEAD(I,LABEL),I=1,2) 09440021 IF (LABEL.EQ.1) PRINT 15,NAME1 09450021 IF (LABEL.GT.1.AND.NAME2.EQ.BLANK) PRINT 15,NAME1 09460021 IF (LABEL.GT.1.AND.NAME2.NE.BLANK) PRINT 25,NAME1,NAME2 09470021 PRINT 35,XLOWR,XHIGH 09480021 XT = (XHIGH - XLOWR) / 12. 09490021 SCALE = 10. / XT 09500021 DO 100 I=1,13 09510021 YT = ABS (XLOWR + XT * FLOAT (I - 1)) 09520021 IAXIS(I) = INT (YT - FLOAT ((INT (YT) / 1000) * 1000)) 09530021 100 CONTINUE 09540021 PRINT 45,(IAXIS(I),I=1,13) 09550021 PRINT 55 09560021 IL = LINE - 1 09570021 XHI = XHIGH 09580021 XLO = XLOWR 09590021 RETURN 09600000 20 CHAR(JT) = BLA 09610002 C 09620001 IL = IL + 1 09630003 IT = IL - (IL / 5) * 5 09640004 JT = INT ((VALUE - XLO) * SCALE + 1.5) 09650005 XT = AST 09660006 IF (JT.LT.1.OR.JT.GT.121) XT = HEX 09670007 IF (JT.LT.1) JT = 1 09680008 IF (JT.GT.121) JT = 121 09690009 CHAR(JT) = XT 09700010 IF (IT.EQ.0) PRINT 65,IL,(CHAR(I),I=1,121),IL 09710011 IF (IT.GT.0) PRINT 75,(CHAR(I),I=1,121) 09720012 RETURN 09730013 C 09740014 30 PRINT 55 09750015 PRINT 45,(IAXIS(I),I=1,13) 09760016 RETURN 09770017 5 FORMAT ('1',/,58X,' P L O T S E T ',I2,/,58X, 09780019 $ '--------------------',//,64X,2A4,/) 09790020 C 09800018 15 FORMAT (62X,8HSTATION ,2A4) 09810021 25 FORMAT (54X,8HSTATION ,2A4,12H TO STATION ,2A4) 09820021 35 FORMAT (/,54X,F11.2,6H TO ,F9.2,///) 09830021 45 FORMAT (6X,12(I3,7X),I3) 09840021 55 FORMAT (7X,25(5H› )) 09850021 65 FORMAT (3X,I3,1H>,121A1,1H<,I3) 09860021 75 FORMAT (7X,121A1) 09870021 END 09880021 SUBROUTINE SUMMAR ( MXDM, MXPS, PX, X, VPV, FDI, GEOCO, U, 09890021 1 PD, NTEST, ANORM, NAME, DAR, VAR, AROT, COUNT, SHIFT,NGR,NDLY, 09900021 1 NDR,LOK,JJ) 09910021 C *********************************************************** 09920025 C *********************************************************** 09930025 C ** ** 09940025 C ** THIS FORTRAN ROUTINE WAS MODIFIED BY J.R.ADAMS OF ** 09950025 C ** MCELHANNEY SURVEYING AND ENGINEERING LTD. TO RUN ** 09960025 C ** ON AN IBM COMPUTER USING FORTRAN G1. THE WORK WAS ** 09970025 C ** DONE FOR THE INDONESIAN GOVERNMENT DURING 1980 AND ** 09980025 C ** 1981. JOB # 083506 ** 09990025 C ** ** 10000025 C *********************************************************** 10010025 C *********************************************************** 10020025 C 10030025 C 10040025 C 10050021 DOUBLE PRECISION ANORM,AROT,COUNT,DNORM,DAR,PX,ROT,U,VAR,Y,Z 10060021 DOUBLE PRECISIONA,A1,CR,DIST,DLAT,DLOG,ESQ,GEOCO,H,PI,S12,X 10070021 DOUBLE PRECISION CMIN , CXYZ ,RHO , WE , RADG C 10090021 DIMENSION X(6,MXPS), GEOCO(3,MXPS), NTEST( 1), PX(MXDM,MXDM), 10100021 1ANORM(MXDM,MXDM),NAME(2,1),PD(1), 10110021 2FDI(5,MXPS),DXS(3) , NSUM(4) , LOK(1) 10120021 C 10130021 DIMENSION COUNT(1), CR(3), LR(3), DAR(3,3), VAR(3,3), AROT(3, 3),10140021 1DNORM(3), NGR(1), ROT(3,3), SHIFT(3,MXPS), U(1),NDLY(1), NDR(4) 10150021 DIMENSION DUMMY(10) LOGICAL*1 DATM(12) QUEK REAL LO,LOK 10170021 C 10180021 C COMMON BLOKS 10190021 C 10200021 COMMON /OPT/ ASG,HOR, SIGA, DXS,TE,XP99,RT,ALG,ACR,OUT,RHO, 10210021 1 FO, LO, WE, RADG, SGMT, ITTG, ITRIM, IWRT,IOUT,INPT, NROW 10220021 COMMON /DBL/ A1,ESQ,PI,CMIN ,CXYZ 10230021 C 10240021 COMMON /COUNT/ NGRPTS, MACOL, MCCOL, MROW, IPASSC, IACCC, 10250021 1 IPASS, NREJT, NDFRJ , IELEV ,IVRJ , MDFM ,NPTS , ISEQ , IACC 10260021 C 10270021 DIST(X,Y,Z)= DSQRT(X**2+Y**2+Z**2) 10280021 DABS(XX)=ABS(XX) 10290021 DSIN(XX)=SIN(XX) 10300021 DCOS(XX)=COS(XX) 10310021 C 10320021 C 10330021 CALL GETDAY (DATM) 10340021 65 PRINT 122, DATM 10350021 C 10360021 JROW=0 10370021 KK=NPTS 10380021 ND=MXDM 10390021 NREJT= NREJT + NDFRJ 10400021 IACC=IPASS-NREJT 10410021 WRITE (IOUT,130) IPASS,NREJT,IACC,NGRPTS,KK 10420021 NREJT= NREJT -NDFRJ 10430021 WRITE(IOUT,164) NREJT,NDFRJ,TE, IELEV,IVRJ 10440021 C 10450021 C RESETTING MATRIX ANORM 10460021 C 10470021 SIGM=1. 10480021 IF(MDFM.GT.0.AND.VPV.GT.0.) SIGM= VPV/MDFM 10490021 IF(ITRIM.EQ.10) SIGM= SIGA 10500021 SGM= SQRT(SIGM) 10510021 PRINT 93, VPV, MDFM, SGM 10520021 93 FORMAT(1H0,34H THE SQUARE SUM OF RESIDUALS VPV=,F15.2,/ 10530021 135H ACCUMULATED DEGREES OF FREEDOM =,5X,I10,/ 10540021 248H THE ESTIMATED ST. DEVIATION OF UNIT WEIGHT SO=, F15.3) 10550021 CALL TRINN( PX, MACOL, 2, MXDM) 10560021 DO 256 I=1,MACOL 10570021 DO 256 J=I,MACOL 10580021 256 PX(J,I)= PX(I,J) 10590021 DO 255 I=1,MACOL 10600021 DO 255 J=1,MACOL 10610021 255 ANORM(I,J)= PX(I,J)* SIGM 10620021 C 10630021 C CENTER ANT. POSITION TO CENTRE 10640021 IR= 0 10650021 CALL CENTER(X, GEOCO, A1, ESQ, NGRPTS, SHIFT, -1.) 10660021 WRITE(IOUT,87) 10670021 DO 79 IN=1,NGRPTS 10680021 IF(NTEST(IN).EQ.10) GO TO 79 10690021 IR= IR +1 10700021 NDLY(IR)= IN 10710021 CALL CARGEO(X(1,IN),X(2,IN), X(3,IN), DLAT,DLOG, H, RAN, A1, ESQ) 10720021 GEOCO(1,IN) = DLAT 10730021 GEOCO(2,IN) = DLOG 10740021 GEOCO(3,IN) = H 10750021 DLOG=DLOG/RADG 10760021 DLAT=DLAT/RADG 10770021 SGN=DLAT/DABS(SNGL(DLAT)) M1=DLAT 10790021 N1=(DLAT-M1)*60.*SGN 10800000 S1=((DLAT-M1)*60.*SGN-N1)*60. 10810001 SGN=DLOG/DABS(SNGL(DLOG)) M2=DLOG 10830003 N2=(DLOG-M2)*60.*SGN 10840004 S2=((DLOG-M2)*60.*SGN-N2)*60. 10850005 WRITE(IOUT,85) (NAME(I,IN),I=1,2), NTEST(IN),PD(IN), M1,N1,S1, 10860006 1 M2,N2,S2, H, (X(I,IN),I=1,3) 10870007 CR(1)=PI/2.-DLAT*RADG 10880008 CR(2)=PI-DLOG*RADG 10890009 LR(1)=2 10900010 LR(2)=3 10910011 C ROTATION MATRIX 10920012 CALL DROTAT (CR,LR,2,ROT) 10930013 JK=3 10940014 JL=JK 10950015 DO 67 I=1,JK 10960016 DO 67 J=1,JK 10970017 DAR(I,J)=0. 10980018 VAR(I,J)=0. 10990019 67 AROT(I,J)=0. 11000020 DO 70 I=1,JK 11010021 DO 69 J=1,JK 11020021 69 VAR(I,J)= ANORM(JROW+I,JROW+J) 11030021 70 CONTINUE 11040021 CALL MULTP (JL,JL,JL,3,3,VAR,ROT,DAR) 11050021 CALL MULTP (JL,JL,JL,3,3,DAR,ROT,VAR) 11060021 IDFM=MDFM 11070021 72 CONTINUE 11080021 DO 73 I=1,JL 11090021 STDEV=DSQRT(VAR(I,I)) 11100021 73 WRITE(IOUT,131) (VAR(I,J),J=1,3), STDEV 11110021 78 JROW=JROW+JK 11120021 79 CONTINUE 11130021 IF(IWRT-9) 155,155,156 11140021 155 WRITE(9) NPTS, MACOL, (NDLY(I),I=1,NPTS),((NAME(I,NDLY(J)),I=1,2),11150021 1(GEOCO(I,NDLY(J)), X(I,NDLY(J)),I=1,3),J=1,NPTS) ,((PX(I,J),I= 11160021 2J,MACOL),J=1,MACOL),MDFM, VPV, IPASS, IELEV, NDFRJ,NREJT, IVRJ 11170021 ENDFILE 9 11180021 156 CONTINUE 11190021 SIGM = SQRT(SIGM) 11200021 WRITE(IOUT,165) 11210021 165 FORMAT('0'//' VARIANCE COVARIANCE MATRIX OF X Y Z ( METERS SQ"D11220021 1)'//) 11230021 DO 166 I=1,MACOL 11240021 166 WRITE(IOUT,125) (ANORM(I,J),J=I,MACOL) 11250021 WRITE(IOUT,173) 11260021 173 FORMAT('0',//' CORRELATION MATRIX OF X,Y,Z '//) 11270021 C 11280021 C COMPUTE CORRELATION MATRIX OF X,Y,Z 11290021 C 11300021 DO 167 I=1,MACOL 11310021 A=DSQRT(ANORM(I,I)) 11320021 DO 168 J=I,MACOL 11330021 168 COUNT(J)= ANORM(I,J)/A/DSQRT(ANORM(J,J)) 11340021 167 WRITE(IOUT,126) (COUNT(J),J=I,MACOL) 11350021 126 FORMAT(' ',8(1X,3F5.2)) 11360021 169 FORMAT('0'//' WEIGHT COEFICIENT MATRIX '//) 11370021 WRITE(IOUT,169) 11380021 DO 170 I=1,MACOL 11390021 170 WRITE(IOUT,125) (PX(I,J),J=I,MACOL) 11400021 IF(IWRT.EQ.10) GO TO 181 11410021 IF(ITRIM.GT.5) GO TO 181 11420021 C REWIND SCRATCH TAPE AND OUTPUT INTERMEADIATE RESULTS 11430021 C JMAX= MAX. OF PASSES TO BE PLOTED 11440021 JMAX= ((JJ/5)/10)*10 11450021 IACTIV = 0 11460021 DO 180 JK=1,NGRPTS 11470021 IF(NTEST(JK).EQ.10) GO TO 180 11480021 IACTIV = IACTIV + 1 11490021 JST=1 11500021 S12= DIST( X(4,JK), X(5,JK), X(6,JK)) 11510021 IS12= (S12/1000) 11520021 IS12= IS12*1000 11530021 DO 172 I=1,3 11540021 LR(I)= (X(I+3,JK)/1000) 11550021 172 LR(I)= LR(I)*1000 11560021 DO 20 I=1,4 11570021 20 NSUM(I)=0 11580021 C IPSS= ACCEPTED PASS NO FOR CURRENT SATION 11590021 IPSS=1 11600021 REWIND 17 1161021 183 FORMAT(' ', 5I4,I3,F7.2,4I3,I4,3X,3(F9.1,F6.1),3X,F9.1,F6.1,3X, 11620021 1 F6.2,I7) 11630021 186 FORMAT('1',//' ',18H SUMMARY AT STA ,2A4,72X,12A1//' ', 11640021 147H ACC PASS SAT DAY LOCKON FREQ. SW SE NW NE ELV,8X,1HX,4X,5HST.11650021 2D.,5X,10HY ST.D.,5X,10HZ ST.D.,6X,12HRADIUS ST.D.,5X,2HSO,4X11660021 3,2HDF/) 11670021 GO TO 191 11680021 190 CONTINUE 11690021 JST= J 11700021 191 WRITE(IOUT,186) NAME(1,JK),NAME(2,JK),(DATM(KX),KX=1,12) 11710021 WRITE(IOUT,30) LR,IS12 11720021 30 FORMAT(' ',49X,3(I8,7X),4X,I7) 11730021 ILMT=IACC 11740021 DO 182 J=JST, IACC 11750021 IPG= IPSS- (IPSS/50)*50 11760021 DO 80 I=1,MACOL 11770021 80 U(I)=0. 11780021 KT=MACOL 11790021 READ(17,END=184) 1180021 + IPASS,IST,(LOK(K), (X(KX,K),KX=1,3),DUMMY(K+NGRPTS),NGR(K)11810000 1,NDR(K),DUMMY(K),K=1,NGRPTS),((PX(I,K),K=I,KT),I=1,KT),SIGM,NDF 11820001 192 CONTINUE 11830003 IF(NGR(JK).EQ.0) GO TO 182 11840004 DO 81 I=1,MACOL 11850005 DO 81 K=I,MACOL 11860006 81 PX(K,I)= PX(I,K) 11870007 FOFF1= DUMMY(JK+ NGRPTS) 11880008 IEL= DUMMY(JK)+.5 11890009 LLOK= LOK(JK) 11900010 CALL DIME(LLOK,IDY,IHJ,IMJ) 11910011 S12= DIST(X(1,JK),X(2,JK),X(3,JK)) -IS12 11920012 JT= IACTIV * 3 -2 11930013 CALL MULTP( 3, 3, 1, ND, ND, PX(JT,JT), X(1,JK),COUNT) 11940014 CALL MULTP( 3, 1, 1, ND, ND, COUNT , X(1,JK), DNORM ) 11950015 SS12= DSQRT(DNORM(1))/IS12 II= NDR(JK) 11970017 DO 10 I=1,4 11980018 10 NDR(I)=0 11990019 NSUM(II)= NSUM(II) +1 12000020 NDR(II)= NSUM(II) 12010021 DO 82 I=1,3 12020021 JT=(IACTIV -1)* 3 +I 12030021 IF(IPSS.LE. JMAX) FDI(I+1,IPSS)= X(I,JK) 12040021 COUNT(I)= X(I,JK)- LR(I) 12050021 82 COUNT(I+3)= DSQRT(PX(JT,JT)) 12060021 IF(IPSS.GT.JMAX) GO TO 250 12070021 FDI(1,IPSS)= FOFF1 12080021 FDI(5,IPSS)= DIST(X(1,JK),X(2,JK),X( 3,JK)) 12090021 NN= IPSS 12100021 250 CONTINUE 12110021 WRITE(IOUT,183) IPSS,IPASS, IST, IDY, IHJ, IMJ, FOFF1, NDR,IEL, 12120021 1 (COUNT(I),COUNT(I+3),I=1,3),S12, SS12, SIGM, NDF 12130021 IPSS= IPSS +1 12140021 IF(IPG.EQ.0) GO TO 190 12150021 182 CONTINUE 12160021 C PLOT DF, X, Y, Z, AND R 12170021 C PLOT ITRIM SETS (DF, X, Y, Z,R) 12180021 184 CONTINUE 12190021 IF(ITRIM.LE.0) GO TO 180 12200021 IPSS= IPSS- 1 12210021 IF(IPSS.GT.0)CALL PPLLOT( ITRIM, FDI,NAME(1,JK),NN, FDI(1,IPSS)) 12220021 180 CONTINUE 12230021 181 RETURN 12240021 C 12250021 85 FORMAT('0',2X,2A4,2X, I2,4X,F6.2,2X,2(I3,I3,F7.3,1X),F10.2,5X,3(F112260021 12.2,3X),/) 12270021 125 FORMAT(' ',(3(3X,3E12.5)) ) 12280021 87 FORMAT('0',//,2X,'ST NO', 5X,'CODE',5X,'WEIGHT',7X,'LAT',9X,'LONG'12290021 1,8X,'HEIGHT',16X,'X',15X,'Y',15X,'Z'/) 12300021 122 FORMAT(1H0,47H P H A S E S O L U T I O N - S U M M A R Y,40X,12310021 ,12A1,//) 12320021 130 FORMAT('0',///,' TOTAL NUM. OF PROCESSED PASSES =',I4,/,' NUMBER12330021 1 OF REJECTED PASSES =',I4,/,' NUMBER OF ACCEPTED PASSES 12340021 2 =',I4,/,' NUMBER OF STATIONS INVOLVED =',I4,/,' NUMBER OF U12350021 3NKNOWN STATIONS =',I4,//) 12360021 131 FORMAT(' ',3X,3E15.5,10X,' ST. DEV.=',F10.2) 12370021 164 FORMAT('0' /' REJECTED ON 99 PER CENT =',I4/2X,'REJECTED ON ZER12380021 10 DEG FRDM =',I4/2X,'REJECTED ON ',F4.1,' DEG ELEV =' ,I4//' DOPPL12390021 2ERS REJECTED ON 99 PER CENT =',I4) 12400021 END 12410021 SUBROUTINE RINPT(MXPS,MXDOP,ND, DOP,NDESC,ISHRT,IEND,NOBR, 12420021 1 KDPSI, NG, DMET, LOK, NSEW, MGR, INJ, NTEST,LOKPR) 12430021 C *********************************************************** 12440025 C *********************************************************** 12450025 C ** ** 12460025 C ** THIS FORTRAN ROUTINE WAS MODIFIED BY J.R.ADAMS OF ** 12470025 C ** MCELHANNEY SURVEYING AND ENGINEERING LTD. TO RUN ** 12480025 C ** ON AN IBM COMPUTER USING FORTRAN G1. THE WORK WAS ** 12490025 C ** DONE FOR THE INDONESIAN GOVERNMENT DURING 1980 AND ** 12500025 C ** 1981. JOB # 083506 ** 12510025 C ** ** 12520025 C *********************************************************** 12530025 C *********************************************************** 12540025 C 12550025 C 12560025 DOUBLE PRECISION A1,ESQ,PI,DOP,CMIN,CXYZ,ORBIT,DT,ISHRT,FS 12570021 DOUBLE PRECISION RHO,WE,RADG 12580021 C 12590021 C INPUT SUBROUTINE 12600021 DIMENSION NDESC(4), DMET(3,MXPS), CXYZ(8,3), LOK(1),NSEW( 1) 12610021 DIMENSION DOP(MXDOP,MXPS),ISHRT(MXDOP,MXPS), ND ( 1),ORBIT( 7) 12620021 1,NOBR(2,1),KDPSI(MXDOP,MXPS),NG(1),DT(4),DXS(3),MGR(1),IFT(4) 12630021 2 ,NTEST(1) ,N(2) 12640009 C 12650010 REAL LOK 12660011 C COMMON BLOCKS 12670012 C 12680013 COMMON /OPT/ ASG,HOR, SIGA, DXS, TE,XP99,RT,ALG,ACR, 12690014 1OUT, RHO, FO,XLO, WE, RADG, SGMT, ITTG, ITRIM, IWRT,IOUT,INPT,12700015 2 NROW 12710016 COMMON /DBL/ A1,ESQ,PI,CMIN,CXYZ 12720021 C 12730017 COMMON /COUNT/ NGRPTS,MACOL,MCCOL,MROW,IPASSC,IACCC, IPASS, NREJT,12740018 1NDFRJ, IELEV, IVRJ, MDFM, NPTS, ISEQ, IACC 12750019 C 12760020 C 12770021 COMMON /SAT/ ORBIT,FS,DT,LOKR,IMIN,IST, NSTA,NDP,BIT,IT2,CPA 12780021 DATA IFT/4H ,3*4HBFIT / 12790021 C 12800021 C 12810021 C NDESC - FIGURE NAME AND SOURCE OF ORBIT, IST- SAT, LOKR- REF.TIME,12820021 C IMIN- TIME SPAN OF ORB.FIT(MIN), INJ-INJECTION ,FS- SAT FQCY, 12830021 7 IACTIV = 0 12840021 IEND=0 12850021 DO 1 J=1,MXDOP 12860021 NG(J)=0 12870021 DO 1 I=1,NGRPTS 12880021 DOP(J,I)=0. 12890021 KDPSI(J,I)=0 12900021 1 CONTINUE 12910021 READ(INPT,END=64) 12920021 + NDESC, IST,INJ,LOKR,IMIN,FS, NSTA,ORBIT,NROW,CXYZ 12930021 1 ,( MGR(I),I=1,NSTA) 12940021 136 CONTINUE 12950021 C SET REFERENCE TIME LOKPR FOR FDI 12960021 IF(IPASS.EQ.IPASSC) LOKPR= LOKR 12970021 IPASS= IPASS+1 12980021 IF(IACTIV .EQ.0) GO TO 9 12990021 IF((IPASS-IPASSC).GT.20.AND.(IACC-IACCC).LE.2) CALL EXIT 13000021 9 CONTINUE 13010021 CALL DIME(INJ,IDYJ,IHRJ,IMNJ) 13020021 CALL DIME(LOKR,IDYR,IHR,IMNR) 13030021 IF(ITTG.NE.5) WRITE(IOUT,91) 13040021 91 FORMAT('1') 13050021 C NROW TROUBLE FLAG .NE.0 BAD FIT(B-FIT) FOR VAR PARAMETERS 13060021 I= NROW +1 13070021 WRITE(IOUT,92) IST, IDYR,IHR,IMNR,NDESC, IMIN, IHRJ, IMNJ, FS, 13080021 1IFT(I), IPASS 13090021 92 FORMAT(1H0,//,10X,I4,I7,I5,I3,5X,7H FIG ^,2A4,4X,5HORB^ ,2A4, 13100021 12X,5HSPAN;,I3,2X,7HINJECT;,2I3,2X,7HS-FQCY;,F8.1,3X,A5,4X, 5HPASS=13110021 2,I4) 13120021 C EXTENDED OUTPUT FORMAT 13130021 IF(ITTG.NE.1) GO TO 3 13140021 WRITE(IOUT,93) 13150021 WRITE(IOUT,94) CXYZ 13160021 93 FORMAT(1H0,10X,35HCHEBYSHEV COEFFICIENTS OF X,Y,Z FIT) 13170021 95 FORMAT(1H0,10X,39HDOPPLER COUNTS AND TIME INTERVALS( MIN)) 13180021 94 FORMAT(1H ,8F13.2) 13190021 3 CONTINUE 13200021 DO 2 I=1,NSTA 13210021 ND(I)=0 13220021 NDD= MGR(I) 13230021 C IF MGR.NE.0 READ MGR DOPPLERS AT THE CURRENT STATION I 13240021 IF(MGR(I).EQ.0) GO TO 2 13250021 READ(INPT,END=64) N,LOK(I), (DMET( K,I), K=1,3),NSEW(I),FDI 13260021 1,IT,( DOP(J,I),J=1,NDD),(ISHRT(J,I),J=1,IT) 13270021 135 CONTINUE 13280021 C CHECK FOR STATION SEQUENCE 13290021 IF(N(1).NE.NOBR(1,I)) GO TO 55 13300022 IF(N(2).EQ.NOBR(2,I)) GO TO 5 13310021 55 WRITE(IOUT,20) N,(NOBR(J,I),J=1,2) 13320021 CALL EXIT 13330021 5 CONTINUE 13340021 20 FORMAT(1H0,31H **** WRONG STATION NAME TAPE5=,2A4,13H COORD CARD =13350021 1,2A4) 13360021 C REINITILIZE ISHRT(1) 13370021 IF(IT.EQ.0) ISHRT(1,I)= DT(1) 13380021 IF(ITTG.NE.1) GO TO 4 13390021 WRITE(IOUT,97) (NOBR(K,I),K=1,2), 13400021 + LOK(I),(DMET(K,I),K=1,3),NSEW(I), FDI 13410021 1,NDD, IT ,ORBIT 13420021 WRITE(IOUT,95) 13430021 WRITE(IOUT,96) (DOP(J,I), ISHRT(J,I),J=1,NDD) 13440021 4 CONTINUE 13450021 97 FORMAT(1H0,2A4,F14.6,F8.1,2F5.1,I3,F8.2,2I3,F10.4,6F10.6) 13460021 96 FORMAT(1H ,(20X,4(F12.1,F12.7))) 13470021 IF( NDD.EQ.0) GO TO 2 13480021 DO 6 J=1,NDD 13490021 IF(DOP(J,I).EQ.0) GO TO 6 13500021 KDPSI(J,I)= J-((J-1)/IT2)*IT2 13510021 ND(I)= ND(I)+1 13520021 NG(J)=NG(J)+1 13530021 6 CONTINUE 13540021 IF(NTEST(I).NE.10) GO TO 8 13550021 ND(I)= 0 13560000 GO TO 2 13570001 8 IACTIV = IACTIV +1 13580002 2 CONTINUE 13590003 IF(IACTIV .EQ.0) GO TO 7 13600004 GO TO 101 13610005 64 IEND=1 13620006 101 RETURN 13630007 END 13640008 SUBROUTINE FORM( MYD, MXDOP2,MXROW, A, C, W, X, 13650009 1HIO, FD , SK, DOP, KDPSI, ISHRT, NTEST, PD, LOKPR,NAME, 13660010 2ELEV, NGR, V, DMET, XS, LOK, NDLY, IC, MGR, NG) 13670011 C *********************************************************** 13680025 C *********************************************************** 13690025 C ** ** 13700025 C ** THIS FORTRAN ROUTINE WAS MODIFIED BY J.R.ADAMS OF ** 13710025 C ** MCELHANNEY SURVEYING AND ENGINEERING LTD. TO RUN ** 13720025 C ** ON AN IBM COMPUTER USING FORTRAN G1. THE WORK WAS ** 13730025 C ** DONE FOR THE INDONESIAN GOVERNMENT DURING 1980 AND ** 13740025 C ** 1981. JOB # 083506 ** 13750025 C ** ** 13760025 C *********************************************************** 13770025 C *********************************************************** 13780025 C 13790025 C 13800025 C 13810012 C SUBROUTINE FORMS DESIGN MATRICES A, C 13820013 C 13830014 DIMENSION A(MXROW,3) ,C(MXROW,MYD), W( 1), SK( 1),DOP( 1),FD(5), 13840015 1 KDPSI( 1), ISHRT( 1), ELEV( 1), NG(1) , X(3), DT(4), CR(4), LR(4)13850016 2, ROT(3,3), AX(3), CXYZ( 8,3), XS(MXDOP2,6), DMET(3),V(1) 13860017 3 , NAME(2),TIIN(4) 13870017 REAL LO, LOK 13880018 C 13890019 C COMMON AREAS 13900020 C 13910021 COMMON /OPT/ ASG, HOR, SIGA, DXS, DYS, DZS, TE, XP99, RT,13920021 1ALG, ACR, OUT, RHO, FO, LO, WE, RADG, SGMT, ITTG, ITRIM, IWRT,13930021 2 IOUT, INPT, NROW 13940021 COMMON /DBL/ A1,ESQ,PI,CMIN,CXYZ 13950021 C 13960021 COMMON /COUNT/ NGRPTS,MACOL,MCCOL,MROW,IPASSC,IACCC, IPASS, NREJT,13970021 1NDFRJ, IELEV, IVRJ, MDFM, NPTS, ISEQ, IACC 13980021 C 13990021 C 14000021 COMMON /SAT/ DTP,AN,EXC, OMEGA, COMEGA,GLAM, DINCL,FS 14010021 + ,DT,LOKR,IMIN,IST, NSTA,NDP, BIT,IT2,CPA 14020021 C 14030021 DOUBLE PRECISION FG,ONE,FDIFR,TX,DUM,T,FMK,EK,DTR 14040021 DOUBLE PRECISION A,C,DIS,DIST,ROT,V,W,SK,CMIN,CXYZ 14050021 DOUBLE PRECISION A1,AX,CR,ESQ,HIO,PI,X,XS,DOP,DT,ISHRT 14060021 + ,DTP,AN,EXC,OMEGA,COMEGA,GLAM,DINCL 14070021 + ,DWE,XLO,FS,RHO,RADG,WE 14080021 C 14090021 DATA LR/3,3,1,3/ 14100021 DATA DWE /4.3752691D-3/ 14110021 C 14120021 DMIN=(LOK- LOKPR )/1440 14130021 C ACTUAL NO. OF DOPPLERS 14140021 MADOP2= MGR +1 14150021 FDI= FD(1) 14160021 FDRIFT= FD(3) 14170021 FDIN= FD(5) 14180021 AIRPR= DMET(1) 14190021 FG = 4.D8 + FDI 14200021 XLO= 2.9979250D8 / FG 14210021 LO = XLO 14220021 TEST= 4.4* FD(2)*60. 14230021 C DISABLE MISCLOSURE TEST WHILE IN SIMULATION MODE 14240021 IF(ITRIM.EQ.10) TEST= 1.E 07 14250021 SWTCH=1.000000E 00 14260021 ONE =(1.D0 - FDI / 4.D8) 14270021 C NTEST= 1 SATELLITE TIME FRAME BR 14280021 C NTEST= 0 RECEIVER TIME FRAME CBR 14290021 IF(NTEST.EQ.0 ) GO TO 13 14300021 ONE = .999998020837D0 14310015 SWTCH= 0. 14320021 13 CONTINUE 14330021 IF( FS.EQ.0) FS= 32000. 14340021 FDIFR = (FDI + FS) * 60.D0 14350021 KKK= KKK+1 14360021 6 CONTINUE 14370021 TT= DMET(2) + 273.13 14380021 EVP= DMET(3) 14390021 IF(AIRPR.EQ.0.) AIRPR=1014./DEXP((HIO)*1.207E-04) 14400021 C DRY COMPONENT HEIGT HD= 40 136 + 148.72(T-273.13) 14410021 DKD= (2.3081E-03 - (7.5793E-03+1.552E-05*HIO )/TT)*AIRPR 14420021 C WET COMPONENT CONSTANT DKW HW=11 000 14430021 DKW= ( 821.2-0.0746*HIO )*EVP /TT/TT 14440021 T= NDLY/.60D8 14450021 IB =IC 14460021 DTX=0. 14470021 P= SQRT(PD) 14480021 DO 23 J=1,MGR 14490021 DTX= (SK(J) -SK(J+1))*SWTCH/CMIN 14500021 TX = ISHRT(J) * ONE + DTX 14510021 DTK =(T - CPA) 14520021 IF(KDPSI(J).EQ.0) GO TO 24 14530021 CCC=0. 14540021 C SAASTAMOINEN TR. REFRACTION ASG=2 MODEL 2 14550021 IF(ASG.EQ.2) CCC= TRRNRC(ELEV(J),ELEV(J+1),AIRPR,EVP,TT, HIO)/LO 14560021 C HOPPFIELD SIMPL. TROP. FORMULA MODEL 1 14570021 IF(ASG.EQ.1) CCC=TRREFR(ELEV(J),ELEV(J+1),DKW,DKD)/LO 14580021 25 IHO= 4 14590021 IHK= 5 14600021 IHD= 6 14610021 C FREQ. OFFSET DERIVATE 14620021 C(IC,IHO) = TX * P 14630021 C(IC,IHD) = CCC/100.* P 14640021 C ABSOLUTE TERM W 14650021 DUM = FDIFR + .83D-3*DTK*FDRIFT 14660021 W(IC) = -(SK(J)-SK(J+1)) /XLO + TX*DUM - DOP(J) + CCC 14670021 C 14680021 TESTD= TEST*TX 14690021 IF(DABS(W(IC)).LT.TESTD ) GO TO 9 14700021 PRINT 28, NAME,J,W(IC) 14710021 28 FORMAT(1H ,33H **** MISCLOSURE REJECTED AT ST.=,2A4,10H DOP NO= ,14720021 1I2,5H W= ,E15.5) 14730021 NGR =NGR-1 14740021 NG(J)=NG(J)-1 14750021 KDPSI(J)=0 14760021 DOP(J)= 0. 14770000 GO TO 24 14780001 C 14790002 C FORMING C MATRIX 14800003 C 14810004 9 CONTINUE 14820005 W(IC)= W(IC)* P 14830006 DIS= SK( J)*XLO/P 14840007 DIST=SK(J+1)*XLO/P 14850008 DO 10 K=1,3 14860009 AX(K)=0. 14870010 C(IC, K )=(X(K) - XS(J,K))/DIS 14880011 C(IC+1, K)=-(X(K)- XS(J+1,K))/DIST 14890012 C FORMING MATRIX A 14900013 10 A(IC,K)=-C(IC, K) - C(IC+1, K) 14910014 C FORM MATRIX C 14920015 C 14930016 C(IC,IHK)=0. 14940017 DO 11 K=1,2 14950018 JK= J-1+K 14960019 DTR= DTP + T + (K-1)*TX 14970020 FMK = AN * DTR 14980021 EK= FMK +EXC*DSIN(FMK) 14990021 CR(1)=-EK 15000021 CR(2)= -OMEGA 15010021 CR(3)= -DINCL 15020021 CR(4)= -COMEGA+GLAM +DWE*DTR 15030021 NO=4 15040021 CALL DROTAT (CR,LR,NO,ROT) 15050021 ICC=IC-1+K 15060021 DO 11 JBC=1,3 15070021 C(IC,IHK)=C(IC,IHK)+XS(JK,JBC+3 )*C(ICC,JBC) /60.D 04 15080021 DO 11 JBB=1,3 15090021 AX(JBC)= AX(JBC) + ROT(JBB,JBC)*C(ICC,JBB) 15100021 11 CONTINUE 15110021 DO 12 JBC=1,3 15120021 12 C(IC,JBC)= AX(JBC) 15130021 IC=IC+1 15140021 24 T= T+ TX 15150021 23 CONTINUE 15160021 IE= IC-1 15170021 C EXTENDED OUTPUT 15180021 IF(ITTG.NE.1) GO TO 18 15190021 WRITE(IOUT,14) 15200021 14 FORMAT(1H0, 10X,69HSATELLITE X,Y,Z AND XDOT,YDOT,ZDOT AT 2MIN INT.15210021 1 STARTING AT LOCK TIME) 15220021 WRITE(IOUT,15)((XS(I,J),J=1,6),I=1,MADOP2) 15230021 15 FORMAT( (5X,3F15.1,5X,3F15.2)) 15240021 WRITE(IOUT,16) 15250021 16 FORMAT(1H0,49H DESIGN MATRICES A, C, W PREMULTIPLIED BY SQRT(P)) 15260021 18 CONTINUE 15270021 NNN=(NGR**2 +NGR )/2 15280021 C 15290021 C NOTE SK AND V ARE ADDRESSES TO BE OVERLAYED IN BLANK COMMON BLOC15300021 C 15310021 IF(RT.EQ.0.) GO TO 19 15320021 CALL WTFORM( MXROW,3,MYD, MGR ,NNN, 3, MYD, A(IB,1),C(IB,1),15330021 1W(IB), SK, KDPSI, 1, 1., RT, V, NGR ) 15340021 19 CONTINUE 15350021 IF(ITTG.EQ.1)WRITE(IOUT,17)((A(I,J),J=1,3),(C(I,J),J=1,MYD),W(I),I15360021 1=IB,IE) 15370021 17 FORMAT(1H ,(5X,3F11.6,5X,6F11.6,5X,F11.2)) 15380021 RETURN 15390021 END 15400021 SUBROUTINE STCORD( ISHRT, MXDOP2, X, NDLY, FDI, SK, 15410021 1 NG,NGR, KDPSI , ELMAX, ELEV , XS,LOK,GCO, DOP, CS,NAME,MGR,NTEST)15420021 C *********************************************************** 15430025 C *********************************************************** 15440025 C ** ** 15450025 C ** THIS FORTRAN ROUTINE WAS MODIFIED BY J.R.ADAMS OF ** 15460025 C ** MCELHANNEY SURVEYING AND ENGINEERING LTD. TO RUN ** 15470025 C ** ON AN IBM COMPUTER USING FORTRAN G1. THE WORK WAS ** 15480025 C ** DONE FOR THE INDONESIAN GOVERNMENT DURING 1980 AND ** 15490025 C ** 1981. JOB # 083506 ** 15500025 C ** ** 15510025 C *********************************************************** 15520025 C *********************************************************** 15530025 C 15540025 C 15550025 C 15560021 C SUBROUTINE COMPUTES SAT. STATE VECTOR XS 15570021 C 15580021 DOUBLE PRECISION Y,TK,CHEB,DCH,DIST,T,DET,ONE 15590021 DOUBLE PRECISION CS,XN,YN,ZN,XU,YU,ZU 15600021 DOUBLE PRECISION A1,ESQ,GCO,PI,X,XS,FS 15610021 DOUBLE PRECISION ALAT,H,SS,DOP,SK,CMIN,CXYZ,ORBIT,DT,ISHRT 15620021 DOUBLE PRECISION RHO,WE,RADG 15630021 C 15640021 DIMENSION ISHRT(1) , SK(1), KDPSI(1), ELEV(1), XS(MXDOP2,6),DT(415650021 1),CS(1),CHEB(8),CXYZ(8,3),DOP(1),DCH( 8), ORBIT( 7) 15660021 DIMENSION NG(1) ,DXS(3), DX(3), GCO(3) , X(3) , NAME(2) 15670021 REAL LO, LOK 15680021 C 15690021 C COMMON AREAS 15700021 C 15710021 COMMON /OPT/ ASG,HOR, SIGA, DXS, TE,XP99,RT,ALG,ACR, OUT, 15720021 1RHO, FO, LO, WE, RADG, SGMT, ITTG, ITRIM, IWRT,IOUT,INPT,NROW 15730021 COMMON /DBL/ A1,ESQ,PI ,CMIN,CXYZ 15740021 C 15750021 COMMON /SAT/ ORBIT,FS,DT,LOKR,IMIN,IST, NSTA,NDP,BIT,IT2,CPA 15760021 C 15770021 ASIN(XX) = ARSIN(XX) 15780021 DIST(XN,YN,ZN,XU,YU,ZU) = DSQRT((XU-XN)**2+(YU-YN)**2+(ZU-ZN)**2) 15790021 C COMPUTE APPROXIMATE COORD. XS , BIT=BIT LENGTH IN MIN 15800021 C CMIN SPEED OF LIGHT IN MIN 15810021 MADOP2= MGR +1 15820021 DO 10 J=1,MADOP2 15830021 10 SK(J)=0. 15840021 ONE = 1.D0 - FDI / 4.D8 15850021 IB=1 15860021 SWTCH=1.000000E 00 15870021 C NTEST= 0 RECEIVER TIME FRAME CBR 15880021 C NTEST= 1 SATELLITE TIME FRAME BR 15890021 IF(NTEST.EQ.0 ) GO TO 14 15900021 C ASSUME SAT TIME FRAME 15910021 ONE = .99999802083D0 15920021 SWTCH= 0. 15930021 14 CONTINUE 15940000 NJET=0 15950001 DET = 0.D0 15960002 5 TK = 2.D0 / (IMIN-1) 15970003 T = DBLE(LOK) - LOKR +NDLY / .6D8 15980004 2 FORMAT(1H ,30X,27H**** WARNING **** DOP NO= ,I3,8H ST NO= ,2A4, 15990005 128H OUTSIDE CHEBYSHEV FIT SPAN) 16000006 C DCC= COS*COS, DSC =COS*SIN DS=SIN OF LATITUDE AND LONGITUDE 16010007 SS= DSQRT(X(1)**2+X(2)**2) 16020008 ALAT= GCO(1) 16030009 H = GCO(3) 16040010 SS= DSQRT(SS**2+(X(3)+(SS/DCOS(ALAT)-H)*ESQ*(DSIN(ALAT)))**2) 16050011 DCC= X(1) /SS 16060012 DCS=X(2) /SS 16070013 DS= SQRT(1.00000E 00- DCC**2- DCS**2)*ALAT/DABS(ALAT) 16080014 DO 1 J=1,MADOP2 16090015 T= T +DET 16100016 Y = T *TK - 1.D0 16110017 CALL CHEBY(CHEB,DCH,1,8,Y,TK) 16120018 DO 3 L= 1,3 16130019 XS(J,L)= CXYZ(1,L)- DXS(L) 16140020 XS(J,L+3) = 0. 16150021 DO 3 K=2,8 16160021 XS(J,L+3) = XS(J,L+3 ) + DCH(K)* CXYZ(K,L) 16170021 3 XS(J,L)= XS(J,L) +CHEB(K)*CXYZ(K,L) 16180021 SK(J)= DIST(X(1),X(2),X(3) ,XS(J,1),XS(J,2),XS(J,3)) 16190021 C 16200021 C CORRECT FOR EARTHS ROTATION DURING RECEIVER PROPAGATION DELAY 16210021 C 16220021 DLY =(NDLY/.6D08 + SK(J)/CMIN) * WE 16230021 XS(J,1) = XS(J,1) + DLY*X(2) 16240021 XS(J,2) = XS(J,2) - DLY*X(1) 16250021 IF(SWTCH.EQ.0.) GO TO 15 16260021 C 16270021 DO 13 K=1,3 16280021 13 XS(J,K)= XS(J,K) + XS(J,K+3 )*(SK( 1)-SK(J))/CMIN 16290021 SK(J)= DIST(X(1),X(2),X(3) ,XS(J,1),XS(J,2),XS(J,3)) 16300021 15 CONTINUE 16310021 DO 16 K=1,3 16320021 16 DX(K)= XS(J,K) -X(K) 16330021 EL = (DX(1)* DCC +DX(2)* DCS +DX(3)* DS )/SK(J) 16340021 CS( J)=(DX(1)*XS(J,4) + DX(2)*XS(J,5 ) +DX(3)* XS(J,6))/SK(J) 16350021 CS(J)= CS(J)/ DIST(0,0,0,XS(J,4),XS(J,5),XS(J,6)) 16360021 C 16370021 C CALCULATE ELEVATIONS AND AND TEST AGAINST HOR 16380021 ELEV(J)= ASIN(EL)/RADG 16390021 C CUT OFF TEST 16400021 IF( J.EQ.1) GO TO 12 16410021 K=J-1 16420021 IF(ELEV(K).GE.HOR.AND.ELEV(J).GT.HOR) GO TO 9 16430021 IF(KDPSI(K).EQ.0) GO TO 9 16440021 KDPSI(K) =0 16450021 DOP(K)= 0. 16460021 NGR= NGR- 1 16470021 NG(K)= NG(K)-1 16480021 NJET=NJET+1 16490021 9 CONTINUE 16500021 C COMPUTE EL MAX 16510021 IF(ELEV(J)- ELEV(K) ) 11,18,18 16520021 11 IF(IB.GE.2) GO TO 12 16530021 IF(J.GT.2) GO TO 17 16540021 ELMAX= ELEV(1) 16550021 CPA= T-DET 16560021 IB=2 16570021 GO TO 12 16580021 17 ELMAX= (ELEV(J) -ELEV(J-2))**2/(8.*(2.*ELEV(K) -ELEV(J)-ELEV(J-2))16590021 1)+ELEV(K) 16600021 CPA=T-CS(J )/(CS(J )-CS(K ))*(DET-(SK(J)-SK(K))/CMIN) 16610021 IB=IB+1 16620021 GO TO 12 16630021 18 ELMAX= ELEV(J) 16640021 12 DET= +ISHRT(J) * ONE 16650021 IF(KDPSI(J).NE.0.AND.T.LT.0.) WRITE(IOUT,2) J, NAME 16660021 1 CONTINUE 16670021 TSPAN= T- IMIN +1 16680021 IF( KDPSI(K).NE.0.AND.TSPAN.GT.0.) WRITE(IOUT,2) K, NAME 16690021 6 CONTINUE 16700021 RETURN 16710021 END 16720021 SUBROUTINE REJECT( V, KDPSI, NG, DOP, MXDOP, IRJ, SGMT, NGR, ISWT,16730021 1IZO,NM) 16740021 C *********************************************************** 16750025 C *********************************************************** 16760025 C ** ** 16770025 C ** THIS FORTRAN ROUTINE WAS MODIFIED BY J.R.ADAMS OF ** 16780025 C ** MCELHANNEY SURVEYING AND ENGINEERING LTD. TO RUN ** 16790025 C ** ON AN IBM COMPUTER USING FORTRAN G1. THE WORK WAS ** 16800025 C ** DONE FOR THE INDONESIAN GOVERNMENT DURING 1980 AND ** 16810025 C ** 1981. JOB # 083506 ** 16820025 C ** ** 16830025 C *********************************************************** 16840025 C *********************************************************** 16850025 C 16860025 C 16870025 DOUBLE PRECISION DOP 16880021 C 16890021 C DATA SCREENING SUBROUTINE 16900021 C 16910021 C V RESIDUALS 16920021 C KDPSI INDEX VECTOR OF USED DOPPLERS 16930021 C NG INDEX VECTOR OF USED DOPPLERS PER INTERVAL MULT STA SOLUT 16940021 C DOP OBSERVATIONS 16950021 C MXDP DIMENSION OF NG, DOP, KDPSI, V PER SATION 16960021 C IVRJ NO OF REJECTIONS 16970021 C SGMT UPPER STAT LIMLIT 16980021 C NGR NO OF OBSERVATIONS PER STAION 16990021 C ISWT SWITCH NE 0 REJECT MAXIMUM RESID ONLY =0 REJECT ALL OUTLAYER17000021 C IZO =1 NO REJECTION ACCURED =0 REJECTION ENCOUNTERED 17010021 C 17020021 DOUBLE PRECISION V 17030021 C 17040021 DIMENSION V(1), KDPSI(1), NG(1), DOP(1), NM(2) 17050021 RMX=0. 17060021 J=1 17070021 DO 1 I=1,MXDOP 17080021 IF(DOP(I).EQ.0.) GO TO 1 17090021 RR = DABS( V(J) ) 17100021 IF(RR.GT.SGMT.AND.ISWT.EQ.0) GO TO 3 17110021 IF(RR.LE.RMX) GO TO 2 17120021 KK=J 17130021 K=I 17140021 RMX=RR 17150021 GO TO 2 17160000 3 DOP(I)=0. 17170001 IZO=0 17180002 KDPSI(I)=0 17190003 NG(I)= NG(I)-1 17200004 NGR= NGR- 1 17210005 IRJ= IRJ +1 17220006 PRINT 5 , V( J),NM, I, IRJ 17230007 2 J= J+1 17240008 1 CONTINUE 17250009 IF(RMX.LE.SGMT.OR. ISWT.EQ.0) GO TO 4 17260010 DOP(K)=0. 17270011 IZO=0 17280012 KDPSI(K)=0 17290013 NG(K)= NG(K)-1 17300014 IRJ= IRJ+1 17310015 NGR= NGR-1 17320016 PRINT 5 , V(KK),NM,K , IRJ 17330017 5 FORMAT(1H ,25H ****REJECTED AT 99: V=,F10.1,8H ST. = ,2A4,8H DO17340018 1P = ,I4, 14H -REJECTION = ,I4) 17350019 4 RETURN 17360020 END 17370021 SUBROUTINE TRINN(A,M,IENTER,MAXP) 17380021 C *********************************************************** 17390025 C *********************************************************** 17400025 C ** ** 17410025 C ** THIS FORTRAN ROUTINE WAS MODIFIED BY J.R.ADAMS OF ** 17420025 C ** MCELHANNEY SURVEYING AND ENGINEERING LTD. TO RUN ** 17430025 C ** ON AN IBM COMPUTER USING FORTRAN G1. THE WORK WAS ** 17440025 C ** DONE FOR THE INDONESIAN GOVERNMENT DURING 1980 AND ** 17450025 C ** 1981. JOB # 083506 ** 17460025 C ** ** 17470025 C *********************************************************** 17480025 C *********************************************************** 17490025 C 17500025 C 17510025 C 13 FEB 78, LAWETZ-CMC 17520021 C FORCE AS MUCH OF THE COMPUTATION AS POSSIBLE TO DP. 17530021 DOUBLE PRECISION FMN,DIV,FMULT,DP,A 17540021 DIMENSION A(MAXP,MAXP) 17550021 IF (M-1) 1,1,2 17560021 1 IF(A(1,1).NE.0.) A(1,1)=DSQRT(1.D0/A(1,1)) 17570021 GO TO 14 17580021 2 DO 6 L=1,M 17590021 FMN=1.D0/A(L,1) 17600021 DIV=A(L,1) 17610021 DO 3 J=2,M 17620021 3 A(L,J-1)=A(L,J)/DIV 17630021 A(L,M)=FMN 17640021 DO 6 J=1,M 17650021 IF (J-L) 4,6,4 17660021 4 FMULT=A(J,1) 17670021 DO 5 K=2,M 17680021 5 A(J,K-1)=A(J,K)-FMULT*A(L,K-1) 17690021 A(J,M)=-FMN*FMULT 17700021 6 CONTINUE 17710021 IF (IENTER-1) 7,7,14 17720021 7 CONTINUE 17730021 A(1,1) = DSQRT(A(1,1)) 17740021 DO 8 I=2,M 17750021 DP=A(1,I) 17760021 8 A(1,I)=DP/A(1,1) 17770021 DO 13 I=2,M 17780021 DIV=0. 17790021 N=I-1 17800021 DO 9 L=1,N 17810021 9 DIV=DIV+A(L,I)*A(L,I) 17820021 A(I,I)=DSQRT(A(I,I)-DIV) 17830021 JL=I+1 17840021 IF (JL-M) 10,10,13 17850021 10 DO 12 J=JL,M 17860021 DIV=0. 17870021 DO 11 L=1,N 17880021 11 DIV=A(L,I)*A(L,J)+DIV 17890021 A(I,J)=(A(I,J)-DIV)/A(I,I) 17900021 12 CONTINUE 17910021 13 CONTINUE 17920021 14 CONTINUE 17930021 RETURN 17940021 END 17950021 SUBROUTINE ADJST( NX, NDPY, NRW, PX, A, C, W, U, ANORM, V, X, PY ,17960021 1Y,VPV, VPV1, PD, PYY, AMTR1,AMTR2,N1, NGR, NBI, NTEST) 17970021 C *********************************************************** 17980025 C *********************************************************** 17990025 C ** ** 18000025 C ** THIS FORTRAN ROUTINE WAS MODIFIED BY J.R.ADAMS OF ** 18010025 C ** MCELHANNEY SURVEYING AND ENGINEERING LTD. TO RUN ** 18020025 C ** ON AN IBM COMPUTER USING FORTRAN G1. THE WORK WAS ** 18030025 C ** DONE FOR THE INDONESIAN GOVERNMENT DURING 1980 AND ** 18040025 C ** ** 18060025 C *********************************************************** 18070025 C *********************************************************** 18080025 C 18090025 C 18100025 C 18110021 DIMENSION PX( NX, NX), C( NRW, NBI), A( NRW, 3), W( NRW), U(NRW),18120021 1 ANORM( NX, NX), X(NX), Y(NDPY), PY (NDPY,NDPY), V(NRW), PD(1) 18130021 DIMENSION PYY(1), AMTR1(NDPY,N1), AMTR2(N1,N1), NGR(1),NTEST(1) 18140021 C SUBROUTINE ADJST 18150021 C PHASE ADJUSTMENT WITH WEIGHTING OF PARAMETERS 18160021 C USAGE 18170021 C PHASE ADJUSTMENT 18180021 C 18190021 DOUBLE PRECISION A,AMTR1,AMTR2,ANORM,C,PX,PY,PYY,RES,U,V,W,Y 18200021 + ,XPX,YPY 18210021 C 18220021 C SUBPROGRAMS REQUIRED 18230021 C MULTP 18240021 C ADMTR 18250021 C REPLC 18260021 C TRIN (GALS SYSTEM SUBR.) 18270021 C 18280021 C COMMON BLOCK AREAS 18290021 C 18300021 COMMON /COUNT/ NGRPTS, MACOL, MCCOL, MROW, IPASSC, IACCC, 18310021 1 IPASS, NREJT, NDFRJ , IELEV ,IVRJ , MDFM ,NPTS , ISEQ , IACC 18320021 C 18330021 ND=N1 18340021 C 18350021 C FORMING MATRIX N OF NORMAL EQNS 18360021 C 18370021 C INITIALIZATION 18380021 JK=0 18390021 C COMPUTE SIGN OF PD 18400021 SIGN= PD(1)/ ABS(PD(1)) 18410021 DO 2 I=1,MCCOL 18420021 DO 1 K=1,MACOL 18430021 1 AMTR2(I,K)=0. 18440021 DO 2 J=1,ND 18450021 AMTR1(I,J)=0. 18460021 2 CONTINUE 18470021 DO 3 I=1,NDPY 18480021 Y(I)=0. 18490000 3 PYY(I)= PY(I,I) 18500001 C 18510002 C START STATION DO-LOOP 18520003 C 18530004 NBO= NBI-3 18540005 IOUT= 3 18550006 IST=1 18560007 IWR=1 18570008 M=0 18580009 DO 8 I=1,NGRPTS 18590010 IF(NTEST(I).NE.10) M= M+1 18600011 IF(NGR(I).EQ.0) GO TO 8 18610012 K= 3*M -2 18620013 L= NGR(I) 18630014 KC=4 +NBO*IST-NBO 18640015 CALL MULTP( L, 3, 3, NRW, NX, A(IWR,1),A(IWR,1),ANORM(K,K)) 18650016 C ANORM =A›A 18660017 C 18670018 C ANORM MATRIX OF NORMAL EQ›NS OF ORBITAL MODE 18680019 C ANORM =PX+A›A 18690020 C FORM U - TEM 18700021 C 18710021 CALL MULTP( L, 3, 1 ,NRW, NX, A(IWR,1), W(IWR), U(K) ) 18720021 CALL MULTP( L,3 ,NBO ,NRW,NDPY,C(IWR, 1), C(IWR, 4),AMTR1(1,KC) )18730021 CALL MULTP( L,NBO,NBO ,NRW,NDPY,C(IWR, 4), C(IWR, 4),AMTR1(KC,KC))18740021 IST= IST +1 18750021 IWR= IWR+ NGR(I) 18760021 8 CONTINUE 18770021 CALL ADMTR( MACOL, MACOL,NX, NX, PX, ANORM, SIGN) 18780021 CALL MULTP(MROW , 3, 3, NRW,NDPY,C, C, AMTR1 ) 18790021 DO 10 I=1,3 18800021 DO 10 J=4,MCCOL 18810021 10 AMTR1(J,I)= AMTR1(I,J) 18820021 C U= A›W 18830021 C AMTR1=C›C 18840021 CALL ADMTR( MCCOL, MCCOL, NDPY, NDPY, AMTR1, PY, 1.) 18850021 CALL TRINN( PY, MCCOL, 2, NDPY) 18860021 DO 26 I=1,MCCOL 18870021 DO 26 J=I,MCCOL 18880021 AMTR1(J,I)=0. 18890021 AMTR1(I,J)=0. 18900021 26 PY(J,I)= PY(I,J) 18910021 C PY(PY+C›C)-1 18920021 IST=1 18930021 IWR=1 18940021 M=0 18950021 DO 29 I=1, NGRPTS 18960021 IF(NTEST(I).NE.10) M= M+1 18970021 IF( NGR(I).EQ.0) GO TO 29 18980021 K= 3*M -2 18990021 L= NGR(I) 19000021 KC= 4+ NBO*IST -NBO 19010021 CALL MULTP( L, 3, 3, NRW,NDPY,C(IWR, 1), A(IWR, 1), AMTR1(1, K) ) 19020021 CALL MULTP( L,NBO,3 , NRW,NDPY,C(IWR, 4), A(IWR, 1), AMTR1(KC,K ))19030021 IST= IST +1 19040021 IWR= IWR+ NGR(I) 19050021 29 CONTINUE 19060021 C 19070021 C AMTR1=C›A 19080021 CALL MULTP( MCCOL, MCCOL, MACOL, NDPY, ND, PY, AMTR1, AMTR2) 19090021 C AMTR2=(PY+C›C)-1C›A 19100021 DO 5 I=1,MACOL 19110021 DO 5 J=1,MACOL 19120021 RES=0. 19130021 DO 6 K=1,MCCOL 19140021 6 RES= RES +AMTR2(K,I)* AMTR1(K,J) 19150021 5 ANORM(I,J)= ANORM(I,J) -SIGN*RES 19160021 C ANORN=A›C( " )-1C›A+A›A+PX 19170021 C 19180021 C ANORM - REDUCED NORMALS 19190021 C 19200021 C V= C" W 19210021 CALL MULTP( MROW, 3, 1, NRW, NRW, C, W, V) 19220021 IWR=1 19230021 IST= 1 19240021 DO 7 K=1,NGRPTS 19250021 IF(NGR(K).EQ. 0) GO TO 7 19260021 J= 4+ IST*NBO- NBO 19270021 CALL MULTP( NGR(K), NBO, 1, NRW,NRW, C(IWR, 4),W(IWR), V(J) ) 19280021 IST=IST +1 19290021 IWR= IWR +NGR(K) 19300021 7 CONTINUE 19310021 DO 27 I=1,MACOL 19320021 DO 27 J=1,MCCOL 19330021 27 U(I)= U(I) -SIGN* AMTR2(J,I)*V(J) 19340021 C U= A"W - A"C( PY + C"C)-1 C"W 19350021 C U - ABSOLUTE TERM 19360021 C SAVE NORMALS 19370021 DO 28 I=1,MACOL 19380021 DO 28 J=1,MACOL 19390021 28 AMTR2(I,J)= ANORM(I,J) 19400021 CALL TRINN(ANORM,MACOL,2,NX) 19410021 C 19420021 C THE WEIGHT COEFICIENT MATRIX (THE INVERSE OF NOR. EQN MATRIX) 19430021 C 19440021 C SOLUTION FOR X- VECTOR 19450021 DO 9 I=1,MACOL 19460021 DO 25 J=I,MACOL 19470021 25 ANORM(J,I)= ANORM(I,J) 19480021 X(I)=0. 19490000 DO 9 J=1,MACOL 19500001 9 X(I)=X(I)-ANORM(I,J)*U(J) 19510002 C X=-ANORM-1 U 19520003 C COMPUTE NUISANCE PARAMETERS !Y 19530004 C 19540005 DO 12 I=1,MCCOL 19550006 DO 11 JA=1,MACOL 19560007 11 V(I)=V(I)+AMTR1(I,JA)*X(JA) 19570008 C V = C›W +C›AX 19580009 12 CONTINUE 19590010 CALL MULTP( MCCOL, MCCOL, 1, NDPY,NDPY, PY, V, Y) 19600011 DO 14 I=1, MCCOL 19610012 14 Y(I)= -Y(I) 19620013 C Y = -(PY+C›C)-1(C›W+C›AX) 19630014 C COMPUTE VECTOR V 19640015 IST=1 19650016 IWR=1 19660017 M=0 19670018 DO 17 K=1,NGRPTS 19680019 IF(NTEST(K).NE.10) M= M+1 19690020 IF(NGR(K).EQ.0) GO TO 17 19700021 IWRE= IWR+ NGR(K)-1 19710021 KK= 3*M - 3 19720021 KKK= + NBO*IST -NBO 19730021 DO 16 I=IWR, IWRE 19740021 V(I)=W(I) 19750021 DO 15 J=1,3 19760021 15 V(I)= V(I)+ A(I,J)*X(J+KK) +C(I,J)*Y(J) 19770021 DO 16 J= 4,NBI 19780021 16 V(I)= V(I)+ C(I,J)* Y(KKK+J ) 19790021 IWR= IWR+ NGR(K) 19800021 IST=IST +1 19810021 17 CONTINUE 19820021 C V = AX+ CY+W 19830021 C 19840021 C COMPUTE AND PRINT THE SQUARE SUM OF RESIDUALS 19850021 C 19860021 DO 18 I=1,MROW 19870021 18 VPV1= VPV1+ SIGN*V(I)*W(I) 19880021 C VPV1=-W›*K 19890021 DO 19 I=1,MROW 19900021 19 VPV= VPV+ SIGN* V(I)**2 19910021 XPX=0. 19920021 DO 21 I=1,MACOL 19930021 W(I)=0. 19940021 DO 20 J=1,MACOL 19950021 20 W(I)=W(I)+X(J)*PX(I,J) 19960021 21 XPX=XPX+W(I)*X(I) 19970021 YPY=0. 19980021 DO 23 I=1,MCCOL 19990021 W(I)=0. 20000021 23 YPY= YPY+ PYY(I)*Y(I)**2 20010021 VPV= VPV+ SIGN*(XPX+YPY) 20020021 C VPV=V›PV+X›PX X 20030021 RETURN 20040021 END 20050021 SUBROUTINE CARGEO (X,Y,Z,ALAT,ALOG,H,RAN,A,ESQ) 20060021 C SUBROUTINE CARGEO 20070025 C 20080025 C DESCRIPTION CONVERTS CARTESIAN COORDINATES TO 20090025 C GEODETIC COORDINATES 20100025 C 20110025 C RESTRICTION ONLY ONE CONVERSION PER CALL 20120025 C 20130025 C X,Y,Z -INPUT CARTESIAN COORDINATES 20140025 C 20150025 C ALAT,ALOG,H - OUTPUT GEODETIC LATITUDE,LONGITUDE (IN RADIANTS) 20160025 C AND GEODETIC HEIGT,(IN METERS),RESPECTIVELY 20170025 C 20180025 C *********************************************************** 20190025 C *********************************************************** 20200025 C ** ** 20210025 C ** THIS FORTRAN ROUTINE WAS MODIFIED BY J.R.ADAMS OF ** 20220025 C ** MCELHANNEY SURVEYING AND ENGINEERING LTD. TO RUN ** 20230025 C ** ON AN IBM COMPUTER USING FORTRAN G1. THE WORK WAS ** 20240025 C ** DONE FOR THE INDONESIAN GOVERNMENT DURING 1980 AND ** 20250025 C ** 1981. JOB # 083506 ** 20260025 C ** ** 20270025 C *********************************************************** 20280025 C *********************************************************** 20290025 C 20300025 C 20310025 C 20320021 C 20330021 DOUBLE PRECISION A,ALAT,ALOG,B,BLAT,ESQ,PI,X,Y,Z,H,D,RON 20340021 C 20350021 PI=3.1415926535D00 20360021 D=Z/(DSQRT(X*X +Y*Y )) 20370021 ALAT=DATAN(D) 20380021 1 BLAT=ALAT 20390021 B=DSIN(ALAT) 20400021 RON=A/ DSQRT(1.000000D 00- ESQ*B*B ) 20410021 RAN = RON 20420021 ALAT=D*(1.D00+ESQ*RON*B/Z) 20430021 ALAT=DATAN(ALAT) 20440021 IF (DABS(ALAT-BLAT).LE.1.D-11) GO TO 2 20450021 GO TO 1 20460021 2 ALOG=DATAN(Y/X) 20470021 IF (ALOG) 3,4,4 20480021 3 IF (Y.LT.0.) ALOG=2.D00*PI+ALOG 20490021 IF (Y.GT.0.) ALOG=PI+ALOG 20500021 GO TO 5 20510021 4 IF (Y.LE.0.) ALOG=PI+ALOG 20520021 5 H= Z/D /DCOS(ALAT)-RON 20530021 RETURN 20540021 END 20550021 SUBROUTINE DROTAT (ANGLE,NAXIS,NUM,ROT) 20560021 C *********************************************************** 20570025 C *********************************************************** 20580025 C ** ** 20590025 C ** THIS FORTRAN ROUTINE WAS MODIFIED BY J.R.ADAMS OF ** 20600025 C ** MCELHANNEY SURVEYING AND ENGINEERING LTD. TO RUN ** 20610025 C ** ON AN IBM COMPUTER USING FORTRAN G1. THE WORK WAS ** 20620025 C ** DONE FOR THE INDONESIAN GOVERNMENT DURING 1980 AND ** 20630025 C ** 1981. JOB # 083506 ** 20640025 C ** ** 20650025 C *********************************************************** 20660025 C *********************************************************** 20670025 C 20680025 C 20690025 C 20700021 DOUBLE PRECISION R1, ROT,ANGLE,R2 20710021 C 20720021 C COMPUTES PRODUCT ROTATION MATRIX IN DOUBLE PRECISION 20730021 C INPUTS 20740021 C ANGLE=SET OF UP TO NINE CONSECUTIVE SINGLE AXIS ROTATIONS 20750021 C NAXIS=SET OF AXES OF ROTATION FOR EACH ANGLE 20760021 C NUM=NUMBER OF ROTATIONS INPUTTED 20770021 C OUTPUT 20780000 C ROT=3X3 PRODUCT ROTATION MATRIX 20790001 DIMENSION ROT(3,3), R1(3,3), R2(3,3), ANGLE(NUM), NAXIS(NUM) 20800002 C 20810005 DO 5 K=1,NUM 20820006 C FIND OTHER TWO AXES N2,N3 20830007 N2=NAXIS(K)+1 20840008 IF (N2.GT.3) N2=N2-3 20850009 N3=N2+1 20860010 IF (N3.GT.3) N3=N3-3 20870011 C SET DIAGONAL ELEMENTS 20880012 R1(NAXIS(K),NAXIS(K))=1. 20890013 R1(N2,N2)=DCOS( ANGLE(K)) 20900014 R1(N3,N3)=R1(N2,N2) 20910015 C SET NONZERO OFF-DIAGONAL ELEMENTS 20920016 R1(N2,N3)=DSIN( ANGLE(K)) 20930017 R1(N3,N2)=-R1(N2,N3) 20940018 C SET ZERO OFF-DIAGONAL ELEMENTS 20950019 R1(NAXIS(K),N2)=0. 20960020 R1(NAXIS(K),N3)=0. 20970021 R1(N2,NAXIS(K))=0. 20980021 R1(N3,NAXIS(K))=0. 20990021 C IF FIRST ROTATION, SET ROT=R1 21000021 IF (K.GT.1) GO TO 2 21010021 DO 1 I=1,3 21020021 DO 1 J=1,3 21030021 1 ROT(I,J)=R1(I,J) 21040021 GO TO 5 21050021 C IF NOT FIRST ROTATION SET R2=R1*ROT AND ROT=R2 21060021 2 DO 3 J=1,3 21070021 DO 3 I=1,3 21080021 R2(I,J)=0. 21090021 DO 3 M=1,3 21100021 3 R2(I,J)=R2(I,J)+R1(I,M)*ROT(M,J) 21110021 DO 4 I=1,3 21120021 DO 4 J=1,3 21130021 4 ROT(I,J)=R2(I,J) 21140021 5 CONTINUE 21150021 RETURN 21160021 END 21170021 SUBROUTINE ADMTR( NR, NC, ND1, ND2, A, B, SGN) 21180021 C *********************************************************** 21190025 C *********************************************************** 21200025 C ** ** 21210025 C ** THIS FORTRAN ROUTINE WAS MODIFIED BY J.R.ADAMS OF ** 21220025 C ** MCELHANNEY SURVEYING AND ENGINEERING LTD. TO RUN ** 21230025 C ** ON AN IBM COMPUTER USING FORTRAN G1. THE WORK WAS ** 21240025 C ** DONE FOR THE INDONESIAN GOVERNMENT DURING 1980 AND ** 21250025 C ** 1981. JOB # 083506 ** 21260025 C ** ** 21270025 C *********************************************************** 21280025 C *********************************************************** 21290025 DOUBLE PRECISION A,B 21300021 C 21310021 C BOTH MATRICES MUST BE COMPATABLE,I.E. THE SAME DIMENSION 21320021 C 21330021 C SGN= SIGN B= A+SNG*B 21340021 C ND2 - ROW DIMENSION OF B 21350021 C ND1 - ROW DIMENSION OF A 21360021 DIMENSION A(1), B(1) 21370021 KA=-ND1 21380021 KB=-ND2 21390021 DO 1 J=1,NC 21400021 KA=KA+ND1 21410021 KB=KB+ND2 21420021 DO 1 I=1,NR 21430021 IA=KA+I 21440021 IB=KB+I 21450021 1 B(IB)= A(IA)+ SGN*B(IB) 21460021 RETURN 21470021 END 21480021 SUBROUTINE MULTP (NRA,NCA,NCB,ND1,ND2,A,B,RES) 21490021 C *********************************************************** 21500025 C *********************************************************** 21510025 C ** ** 21520025 C ** THIS FORTRAN ROUTINE WAS MODIFIED BY J.R.ADAMS OF ** 21530025 C ** MCELHANNEY SURVEYING AND ENGINEERING LTD. TO RUN ** 21540025 C ** ON AN IBM COMPUTER USING FORTRAN G1. THE WORK WAS ** 21550025 C ** DONE FOR THE INDONESIAN GOVERNMENT DURING 1980 AND ** 21560025 C ** 1981. JOB # 083506 ** 21570025 C ** ** 21580025 C *********************************************************** 21590025 C *********************************************************** 21600025 C 21610021 DOUBLE PRECISION A,B,RES 21620021 C SUBROUTINE COMPUTES A PRODUCT OF TWO GENERAL MATRICES A›*B 21630021 C NOTE NO. OF ROWS OF A = NO. OF ROWS OF B 21640021 C 21650021 C ND1 - ROW DIMENSION OF A AND B 21660021 C ND2 - ROW DIMENSION OF RES 21670021 DIMENSION A(1), B(1), RES(1) 21680021 KB=-ND1 21690021 KR=-ND2 21700021 DO 1 J=1,NCB 21710021 KA=-ND1 21720021 KB=KB+ND1 21730021 KR=KR+ND2 21740021 DO 1 I=1,NCA 21750021 KA=KA+ND1 21760021 LR=KR+I 21770021 RES(LR)=0. 21780021 DO 1 II=1,NRA 21790021 LA=KA+II 21800021 LB=KB+II 21810021 1 RES(LR)=RES(LR)+A(LA)*B(LB) 21820021 RETURN 21830000 END 21840001 SUBROUTINE WTFORM(NRW,MXDM,MYDM, MXDOP,MXP,MACL,MCCL, A, C, W,21850021 1 P, KDSPI, NSTA, PW, RAT, PT, NGR) 21860021 C *********************************************************** 21870025 C *********************************************************** 21880025 C ** ** 21890025 C ** THIS FORTRAN ROUTINE WAS MODIFIED BY J.R.ADAMS OF ** 21900025 C ** MCELHANNEY SURVEYING AND ENGINEERING LTD. TO RUN ** 21910025 C ** ON AN IBM COMPUTER USING FORTRAN G1. THE WORK WAS ** 21920025 C ** DONE FOR THE INDONESIAN GOVERNMENT DURING 1980 AND ** 21930025 C ** 1981. JOB # 083506 ** 21940025 C ** ** 21950025 C *********************************************************** 21960025 C *********************************************************** 21970025 C 21980021 DOUBLE PRECISION A,C,PT,W,P 21990021 C 22000021 DIMENSION A(NRW,MXDM),C(NRW,MYDM), W( 1), KDSPI( 1), 22010021 1 P(MXP) , PT(MXDOP,MXDOP), PW(1) ,NGR(1) 22020021 C FOR A.X + C.Y -V =W, 22030021 C WE WISH TO PREMULTIPLY BY SQRT(P) 22040021 C ON RETURN - A = SQRT(P) * A 22050021 C - C = SQRT(P) * C 22060021 C - W = SQRT(P) * W. 22070010 C CAUTION, RESIDUALS WILL NOT BE STANDARDIZED. 22080011 KJ=0 22090012 DO 30 I=1,MXP 22100013 30 P(I)=0. 22110014 DO 10 I=1,NSTA 22120015 DO 31 K=1,MXDOP 22130016 DO 31 J=1,MXDOP 22140017 31 PT(K,J)=0. 22150018 M=0 22160019 DO 20 J=1,MXDOP 22170020 IF(KDSPI( J))20,20,9 22180021 9 M=M+1 22190021 PT(M,M)= 1./PW(I) 22200021 IF((J+1)-MXDOP)7,7,20 22210021 7 IF(KDSPI( J+1))20,20,11 22220021 11 PT(M,M+1)= RAT/ PW(I) 22230021 PT(M+1,M)=PT(M,M+1) 22240021 20 CONTINUE 22250021 CALL TRINN(PT,M,1, MXDOP) 22260021 DO 40 K=1,M 22270021 DO 40 J=K,M 22280021 KJ= KJ+1 22290021 40 P(KJ)= PT(K,J) 22300021 10 CONTINUE 22310021 CALL MATMLP(P, A, W, MACL,NRW,NSTA,NGR) 22320021 CALL MATMLP(P, C, PT, MCCL, NRW,NSTA,NGR) 22330021 RETURN 22340021 END 22350021 SUBROUTINE MATMLP( P, A, C, M, MAXP,NSTA,NGR) 22360021 C *********************************************************** 22370025 C *********************************************************** 22380025 C ** ** 22390025 C ** THIS FORTRAN ROUTINE WAS MODIFIED BY J.R.ADAMS OF ** 22400025 C ** MCELHANNEY SURVEYING AND ENGINEERING LTD. TO RUN ** 22410025 C ** ON AN IBM COMPUTER USING FORTRAN G1. THE WORK WAS ** 22420025 C ** DONE FOR THE INDONESIAN GOVERNMENT DURING 1980 AND ** 22430025 C ** 1981. JOB # 083506 ** 22440025 C ** ** 22450025 C *********************************************************** 22460025 C *********************************************************** 22470025 C 22480021 DOUBLE PRECISION A,C,SUM,P 22490021 C 22500021 DIMENSION A(1), C(1), P(1),NGR(1) 22510021 C 22520021 C CALCULATES A= P*A AND C= P*C 22530021 C WHERE 22540021 C MATRIX C IS N BY 1 22550021 C MATRIX A IS N BY M 22560021 C WHERE N IS SUM NGR(I) 22570021 C MATRIX P IS UPPER TRIANGLE MATRIX 22580021 C 22590021 C MAXP - ROW DIMENSION OF A 22600021 C 22610021 IJ=0 22620021 IJJ=0 22630021 I=0 22640021 DO 5 IF=1,NSTA 22650021 KKK= NGR(IF) 22660021 DO 5 II=1,KKK 22670021 IDA= -MAXP 22680021 I= I+1 22690021 IH= I+KKK-II 22700021 DO 2 K=1,M 22710021 SUM=0. 22720021 IDA= IDA+MAXP 22730021 DO 1 J=I,IH 22740021 IJ= IJJ+J -I+1 22750021 JK= IDA+J 22760021 1 SUM= SUM+ P(IJ)* A(JK) 22770021 IK= IDA+I 22780021 2 A(IK)= SUM 22790021 IJJ= IJJ+KKK-II+1 22800021 5 CONTINUE 22810021 IJ=0 22820021 I=0 22830021 DO 4 IF=1,NSTA 22840021 KKK=NGR(IF) 22850021 DO 4 II=1,KKK 22860021 I= I+1 22870021 IH= I+ KKK-II 22880021 SUM=0. 22890021 DO 3 J=I,IH 22900021 IJ= IJ+1 22910021 3 SUM= SUM + P(IJ)*C(J) 22920021 4 C(I)= SUM 22930021 RETURN 22940021 END 22950021 FUNCTION TRRNRC(ELEV1,ELEV2,AIRPR,ES,T,H) 22960021 DOUBLE PRECISION H 22970021 C AIRPR=PRESSURE IN MILLIBARS. 22980021 PI= 3.1415926E 00 22990021 RAD=57.29578E 00 23000021 E1= ELEV1/RAD 23010021 E2= ELEV2/RAD 23020021 C1=DRHOI(AIRPR,T,ES,E1,PI,H) 23030021 C2=DRHOI(AIRPR,T,ES,E2,PI,H) 23040021 TRRNRC=C2-C1 23050021 RETURN 23060021 END 23070021 FUNCTION DRHOI(P,T,E,EL,PI,H) 23080021 DOUBLE PRECISION H,SZ,TZ,Z 23090021 Z=PI/2.-EL 23100021 TZ = DSIN(Z) / DCOS(Z) 23110021 W=16.*4.848E-06*TZ/T*(P+4810.*E/T) 23120021 Z=Z-W 23130000 SZ = 1.D0 / DCOS(Z) 23140001 TZ = SZ * DSIN(Z) 23150002 CONH = 1.16 + H * (-.15412E-3 + .71633E-8*H) 23160003 DRHOI=.002277*SZ*(P+(1255./T+0.05)*E-CONH*TZ*TZ) 23170004 RETURN 23180005 END 23190006 REAL FUNCTION TRREFR(ELE1, ELE2, DKW, DKD) 23200007 DOUBLE PRECISION RAD 23210021 C 23220021 C 23230021 DABS(XX)=ABS(XX) 23240008 C SUBPROGRAM EVALUATES TROPOSPHERIC REFR. CORRECTION 23250009 C USING HOPFIELD›S SIMPLIFIED FORMULAE 23260010 C ELEV ELEVATION ANGLE DEG 23270011 RAD=57.29577951D0 23280012 ANGL2=SQRT(ELE2**2+6.25)/RAD 23290013 ANGLE=SQRT(ELE1**2+6.25)/RAD 23300014 TRREFR= DKD*( 1./SIN(ANGL2)-1./SIN(ANGLE)) 23310015 C WET COMPONENT 23320016 C CALCULATE WET COMPONENT CONSTANT DKW 23330017 ANGL2=SQRT(ELE2**2+2.25)/RAD 23340018 ANGLE=SQRT(ELE1**2+2.25)/RAD 23350019 2 TRREFR=TRREFR+DKW*(1./SIN(ANGL2)-1/SIN(ANGLE)) 23360020 RETURN 23370021 END 23380021 SUBROUTINE CHEBY (CH,DC,M,N,X,XP) 23390021 C *********************************************************** 23400025 C *********************************************************** 23410025 C ** ** 23420025 C ** THIS FORTRAN ROUTINE WAS MODIFIED BY J.R.ADAMS OF ** 23430025 C ** MCELHANNEY SURVEYING AND ENGINEERING LTD. TO RUN ** 23440025 C ** ON AN IBM COMPUTER USING FORTRAN G1. THE WORK WAS ** 23450025 C ** DONE FOR THE INDONESIAN GOVERNMENT DURING 1980 AND ** 23460025 C ** 1981. JOB # 083506 ** 23470025 C ** ** 23480025 C *********************************************************** 23490025 C *********************************************************** 23500025 C 23510021 C 23520021 C WRITTEN BY P. LAWNIKANIS MARCH 1974. 23530021 C 23540021 C 23550021 C ›CHEBY› COMPUTES N-1ST CHEBYSHEV POLYNOMIALS 23560021 C FOR ARGUEMENT ›X› IN VECTOR ›CH›. 23570021 C DERIVATIVE POLYNOMIALS IN ›DC› FOR 23580021 C ›X› DERIVATIVE ›XP› IF ›J› = 1. 23590021 C 23600021 DOUBLE PRECISION CH,DC,TX,X,XP 23610021 C 23620021 DIMENSION DC(N),CH(N) 23630021 IF (N) 10,10,20 23640021 10 RETURN 23650021 C 23660021 20 CH(1) = 1.D0 23670021 CH(2) = X 23680021 IF (N.LT.3) RETURN 23690021 TX = X + X 23700021 DO 100 I = 3,N 23710021 CH(I) = CH(I-1) * TX - CH(I-2) 23720021 100 CONTINUE 23730021 IF (M.NE.1) RETURN 23740021 C 23750021 DC(1) = 0.D0 23760021 DC(2) = XP 23770021 DC(3) = (TX + TX) * XP 23780021 IF (N.LT.4) RETURN 23790021 DO 200 I = 4,N 23800021 J = I - 1 23810021 K = I - 2 23820021 DC(I) = DC(J) * DFLOAT(J) / DFLOAT(K) * TX 23830021 + -DC(K) * DFLOAT(J) / DFLOAT(K-1) 23840021 200 CONTINUE 23850021 RETURN 23860021 C 23870021 END 23880021 SUBROUTINE OPTION ( IN, IEND,MRCN) 23890021 C *********************************************************** 23900025 C *********************************************************** 23910025 C ** ** 23920025 C ** THIS FORTRAN ROUTINE WAS MODIFIED BY J.R.ADAMS OF ** 23930025 C ** MCELHANNEY SURVEYING AND ENGINEERING LTD. TO RUN ** 23940025 C ** ON AN IBM COMPUTER USING FORTRAN G1. THE WORK WAS ** 23950025 C ** DONE FOR THE INDONESIAN GOVERNMENT DURING 1980 AND ** 23960025 C ** 1981. JOB # 083506 ** 23970025 C ** ** 23980025 C *********************************************************** 23990025 C *********************************************************** 24000025 C 24010021 C READS IN AND WRITES OUT OPTIONS, SET DEFAULT VALUES AND CHECKS FOR 24020021 C EOF ON TAPE99=INPUT 24030021 C IN=1 READ + WRITE 24040021 C IN=0 WRITE OUT ONLY 24050021 C IEND=1 EOF ENCOUNTERED ON TAPE99 24060021 C 24070021 DOUBLE PRECISION AD1, B1 24080021 DOUBLE PRECISION A1,ESQ,PI ,CMIN,CXYZ 24090021 DOUBLE PRECISION RHO,WE,RADG 24100021 C 24110021 DIMENSION DXS(3),DATM(3) 24120021 COMMON /OPT/ ASG,HOR, SIGA, DXS,TE,XP99,RT,ALG,ACR,OUT,RHO, 24130021 1 FO, LO, WE, RADG, SGMT, ITTG, ITRIM, IWRT,IOUT,INPT, NROW 24140021 COMMON /DBL/ A1,ESQ,PI, CMIN,CXYZ 24150021 COMMON /COUNT/ NGRPTS, MACOL, MCCOL, MROW, IPASSC, IACCC, 24160021 1 IPASS, NREJT, NDFRJ , IELEV ,IVRJ , MDFM ,NPTS , ISEQ , IACC 24170021 C 24180021 CCARD READER SELECT CODE 24190021 IREAD=5 24200021 IEND=0 24210021 C O P T I O N C A R D READ IN 24220021 C 24230021 C ASG =1 TROP. REFRACTION HOPPFIELD"S SIMPLIFIED MODEL 24240021 C =2 NRC SASTAMONIAN MODEL 24250021 C .NE. 1 TROP. REF. SUPPRESSED 24260021 C HOR - HORIZON CUT-OFF CRITERION (DEGREES) 24270021 C ITRIM = 1 TRIMMING ALL DOPPLER TO KEEP SIMUL. DOPPLER ONLY 24280021 C FROM ALL STA. 24290021 C SIGA - APRIORI VAR FACTOR 24300021 C A1 - SEMIMAJOR AXIS OF REF ELLIPSOID 24310021 C B1 - SEMIMINOR AXIS OF REF ELLIPSOID 24320021 C DXS,DYS,DZS - DATUM SHIFT IN METERS 24330021 C RT - CORRELATION COEFFICIENT(SERRIALLY CORR. P-MATRIX OF DOPPLERS) 24340021 C ALG - ALONG TRK ORBIT CONSTRAIN METERS 1 SIGMA 24350021 C ACR - ACCRR.TRK ORBIT CONSTRAIN METERS 1 SIGMA 24360021 C OUT - OUT - TRK ORBIT CONSTRAIN METERS 1 SIGMA 24370021 C 24380021 C RT - CORRELATION COEFFICIENT WHEN RT=0 IDENTITY WEIGHTING 24390021 C 24400021 IF(IN.EQ.0) GO TO 19 24410021 READ(IREAD,84,END=16) 24420021 + ITTG,ASG, HOR, ITRIM, SIGA, A1, B1,(DXS(I),I=1,3) , 24430021 1 TE, RT, IWRT, ALG, ACR, OUT, MRCN 24440021 84 FORMAT (I2,F3.0,F5.1,3X,I2,F5.1,7F5.1,3X,I2,3F5.1,4X,I1) 134 CONTINUE 24460021 C XP 99 P.C. NOR DIST. LIMITS + APRIORI VAR FOR VEIGHTING 24470021 C CHANGE IF REQUIRED 24480000 XP99=2.576 24490001 IF(ALG.EQ.0.) ALG=26. 24500002 IF(ACR.EQ.0.) ACR= 5. 24510003 IF(OUT.EQ.0.) OUT= 10. 24520004 IF(SIGA.EQ.0) SIGA=1.0 24530005 NDFRJ=0 24540006 IELEV=0 24550007 IF(A1.EQ.0.) A1=135.0 24560008 IF(B1.EQ.0.) B1=750.52 24570009 AD1 = 6378000.D0 + A1 24580021 A1= 6378000. +A1 24590010 B1= 6356000. +B1 24600011 IF(ASG.GT.2.) ASG=0. 24610012 XP99= XP99*SQRT(SIGA) 24620013 IF(TE.EQ.0.) TE=14.5 24630014 181 FORMAT('0',14X,'TE =',F5.1,' PASS ELEV CUT OFF'//15X, 'IWRT =24640015 1',I5,' =1 READ IN (FROM INPUT) PHASE WEIGHT MATRIX +WRITE TAPE 9 24650016 2CONT FILE'/28X,'=9 READ TAPE9 CONT FILE,REWOUND'/28X,'=8 READ TAPE24660017 38 (NOT REWOUND)'/28X,'=10 OUTPUT CLEANED DATA ON TAPE 7'/29X,' 24670018 4IF .GT.9 NO TAPE9 CONT FILE WRITTEN') 24680019 182 FORMAT('0',14X,'RT =',F5.2,' CORRELATION COEFFICIENT'/28X, 24690020 1'=0 IDENTITY MATRIX WEIGHTING') 24700021 AD1 = B1 / AD1 24710021 ESQ = (1.D0 - AD1) * (1.D0 + AD1) 24720021 NREJT=0 24730021 ISEQ=0 24740021 IVRJ=0 24750021 19 CONTINUE 24760021 CALL GETDAY(DATM) 24770021 C PRINT OUT OPTIONS 24780021 PRINT 163 ,DATM 24790021 PRINT 85,ITTG 24800021 PRINT 81, ASG,HOR,ITRIM,SIGA 24810021 PRINT 181,TE, IWRT 24820021 PRINT 182,RT 24830021 WRITE(IOUT,83) ALG,ACR,OUT,A1,B1 24840021 PRINT 82, DXS 24850021 C 24860021 GO TO 17 24870021 16 IEND=1 24880021 17 RETURN 24890021 81 FORMAT(//,15X,'ASG =',F5.0,' (ASG=1 HOPFIELD TR REFRACTION MOD24900021 1EL' /32X,'=2 SAASTAMOINEN MODEL'// 15X,'HOR =',F5.1,' DEGRE24910021 2ES HORIZON CUT-OFF CRITERION'//15X,'ITRIM=',I5,' NO OF SETS OF24920021 3(FDI,X,Y,Z,R) TO BE PLOTTED ,IF.GT.5 NO STATION SUMMARY'/ 24930021 4 28X,'=10 -SIMULATION MODE'//15X, 'SIGM =',F5.1,' APRIORI 24940021 5VARIANCE FACTOR (FOR STAT TESTING)') 24950021 85 FORMAT('0',///,10X,'OPTIONS USED FOR THIS RUN ARE',///,15X,'ITTG =24960021 1',I5, ' =1 EXTENDED OUTPUT',/28X,'=5 SHORT OUTPUT'/28X,'=0 NORMAL24970021 2 OUTPUT') 24980021 82 FORMAT('0',/,10X,'DATUM SHIFT (METERS)',//,15X,'DX=',F8.2,//,15X,'24990021 1DY=',F8.2,//,15X,'DZ=',F8.2) 25000021 83 FORMAT('0',14X,'ALONG=',F5.1,' ORBIT CONSTRAINS 1 SIG.',/,15X,'AC25010021 1CRS =',F5.1,/,15X,'OUT =',F5.1,' METERS', 25020021 2 //,10X,'REFERENCE ELLIPSOID',//,15X,'SEM.MAJ. A=',F10.2,25030021 1 ' M ',15X,' SEM MIN. B=',F10.2,//) 25040021 163 FORMAT('1',' GEODETIC SURVEY OF CANADA PROGRAM G E O D O P (VER25050021 1SION III,OCT.74)',40X, 3A4) 25060021 END 25070021 SUBROUTINE DIME(KTC, IDY,IHR,MIN) 25080021 C *********************************************************** 25090025 C *********************************************************** 25100025 C ** ** 25110025 C ** THIS FORTRAN ROUTINE WAS MODIFIED BY J.R.ADAMS OF ** 25120025 C ** MCELHANNEY SURVEYING AND ENGINEERING LTD. TO RUN ** 25130025 C ** ON AN IBM COMPUTER USING FORTRAN G1. THE WORK WAS ** 25140025 C ** DONE FOR THE INDONESIAN GOVERNMENT DURING 1980 AND ** 25150025 C ** 1981. JOB # 083506 ** 25160025 C ** ** 25170025 C *********************************************************** 25180025 C *********************************************************** 25190025 C KTC TIM IN MIN 25200021 C IDY DAY NO 25210021 C IHR HOURS 25220021 C MIN MINUTES 25230021 IDY= KTC/ 1440 +.5 25240021 IHR=(KTC-IDY*1440)/60 +.5 25250021 MIN=KTC-IDY*1440- IHR*60+.5 25260021 RETURN 25270021 END 25280021 SUBROUTINE GETDAY(DAYT) C LOGICAL*1 DATE(18),DAYT(12) C C CALL SUBROUTINE GDATE IN THE LIBRARY CALL GDATE(DATE) C SELECT DESIRED ELEMENTS DAYT(1) = DATE(6) DAYT(2) = DATE(7) DAYT(3) = DATE(8) DO 10 I = 10,18 L = I-6 10 DAYT(L) = DATE(I) RETURN END