c c alle Subroutinen, die fuer PPM-Ausgleichung gebraucht werden c inclusive IMSL Subroutinen c SUBROUTINE TPMPM (XA,XD,XMUE,XMUES,DMUFOR,T0,T1,YA,YD,YMUE,YMUES) C C EPOCHENUEBERTRAGUNG VON OERTERN UND EIGENBEWEGUNGEN C C EINGABE: XA = ALPHA FUER AUSGANGSEPOCHE C XD = DELTA FUER AUSGANGSEPOCHE C XMUE = E.B. IN ALPHA ZUR AUSGANGSEPOCHE C XMUES = E.B. IN DELTA ZUR AUSGANGSEPOCHE C T0 = AUSGANGSEPOCHE C T1 = ENDEPOCHE C DMUFOR = FORESHORTENING-EFFEKT C AUSGABE: YA = ALPHA FUER ENDEPOCHE C YD = DELTA FUER ENDEPOCHE C YMUE = E.B. IN ALPHA ZUR ENDEPOCHE C YMUES = E.B. IN DELTA ZUR ENDEPOCHE C ALLE WINKEL IM BOGENMASS C ZEITEINHEIT DER EIGENBEWEGUNGEN: JULIANISCHES JAHRHUNDERT C T0, T1 IN JULIANISCHEM DATUM C C IMPLICIT REAL*8 (A-H,O-Z) REAL*8 XM(3,1),ABLA(3,1),ABLD(3,1) C C DT = (T1-T0)/36525.D0 C SD = DSIN(XD) CD = DCOS(XD) SA = DSIN(XA) CA = DCOS(XA) C XMUE0 = DSQRT(XMUE*XMUE*CD*CD + XMUES*XMUES) XMUTOT = XMUE0 + 0.5D0*DMUFOR*DT C C FALLS DER BETRAG DER GESAMTEN EIGENBEWEGUNG EXTREM KLEIN IST, C WIRD SIE ZU NULL GESETZT UND DAS PROGRAM OHNEWEITERE RECHNUNG C SOFORT VERLASSEN C IF(XMUE0.LT.1.D-20) *THEN YMUE =XMUE YMUES=XMUES YA=XA YD=XD RETURN END IF C C BEGINN DER EIGENTLICHEN RECHNUNG C SPSI = (XMUE*CD)/XMUE0 CPSI = XMUES/XMUE0 SMT = DSIN(XMUTOT*DT) CMT = DCOS(XMUTOT*DT) C C BERECHNUNG DES VEKTORS "M-PUNKT MAL X" (SIEHE UNTERLAGEN) C C 'ABLEITUNG DER ORTS-AENDERUNG BERECHNEN', SIEHE UNTERLAGEN C XMUPKT = XMUTOT + 0.5D0*DMUFOR*DT C XM(1,1) = -XMUPKT*(+SD*CA*CPSI*CMT + SA*SPSI*CMT + CD*CA*SMT) XM(2,1) = -XMUPKT*(+SD*SA*CPSI*CMT - CA*SPSI*CMT + CD*SA*SMT) XM(3,1) = -XMUPKT*(-CD*CPSI*CMT + SD*SMT) C C BERECHNUNG DES ORTES FUER DIE ENDEPOCHE C CALL TPMPOS (XA,XD,XMUE,XMUES,DMUFOR,T0,T1,YA,YD) C C BERECHNUNG DES VEKTORS DER ABGELEITETEN RICHTUNGSKOSINUSSE C FUER DIE OERTER ZUR ENDEPOCHE C CALL DLDA (YA,YD,ABLA) CALL DLDD (YA,YD,ABLD) C C TRANSFORMIERTE EIGENBEWEGUNGEN BERECHNEN C CALL SCLPRD (XM,ABLA,YMUE ,3) CALL SCLPRD (XM,ABLD,YMUES,3) YMUE = YMUE/(DCOS(YD)*DCOS(YD)) RETURN END C SUBROUTINE TPMPOS (XA,XD,XMUE,XMUES,DMUFOR,T0,T1,YA,YD) C C EIGENBEWEGUNGSUEBERTRAGUNG VON MITTLEREN OERTERN C C EINGABE: XA = ALPHA ZUR EPOCHE T0 C XD = DELTA ZUR EPOCHE T0 C XMUE = E.B. IN ALPHA ZUR EPOCHE T0 C XMUES = E.B. IN DELTA ZUR EPOCHE T0 C DMUFOR = FORESHORTENING-EFFEKT C T0 = AUSGANGSEPOCHE C T1 = ENDEPOCHE C AUSGABE: YA = ALPHA ZUR ENDEPOCHE C YD = DELTA ZUR ENDEPOCHE C C ALLE WINKEL IM BOGENMASS C ZEITEINHEIT JULIANISCHES JAHRHUNDERT C C IMPLICIT REAL*8 (A-H,O-Z) REAL*8 XM(3,1) C PI = 3.141592653589793D0 C VORBEREITUNGEN DT = (T1-T0)/36525.D0 C SD = DSIN(XD) CD = DCOS(XD) SA = DSIN(XA) CA = DCOS(XA) C XMUE0 = DSQRT(XMUE*XMUE*CD*CD + XMUES*XMUES) XMUTOT = XMUE0 + 0.5D0*DMUFOR*DT C FALLS DER BETRAG DER GESAMTEN EIGENBEWEGUNG EXTREM KLEIN IST, C WIRD SIE ZU NULL GESETZT UND DAS PROGRAM OHNEWEITERE RECHNUNG C SOFORT VERLASSEN C IF(XMUE0.LT.1.D-20) *THEN YA=XA YD=XD RETURN END IF C BEGINN DER EIGENTLICHEN RECHNUNG C C SPSI = (XMUE*CD)/XMUE0 CPSI = XMUES/XMUE0 C SMT = DSIN(XMUTOT*DT) CMT = DCOS(XMUTOT*DT) C C BERECHNUNG DES VEKTOR "M.X" NACH MUELLER, S.115, FORMEL 4.94 C XM(1,1) = -SD*CA*CPSI*SMT - SA*SPSI*SMT + CD*CA*CMT XM(2,1) = -SD*SA*CPSI*SMT + CA*SPSI*SMT + CD*SA*CMT XM(3,1) = +CD*CPSI*SMT + SD*CMT C C BERECHNUNG DER WINKEL AUS DEN RICHTUNGSKOSINUSSEN C CALL ANGLE (XM,YA,YD) RETURN END C SUBROUTINE DLDD (A,D,V) C C BERECHNUNG DES NACH DELTA ABGELEITETEN VEKTORS DER C RICHTUNGSKOSINUSSE C C EINGABE: A, D = ALPHA,DELTA (IM BOGENMASS) C AUSGABE: V = VEKTOR DER RICHTUNGSKOSINUSSE C IMPLICIT REAL*8 (A-H,O-Z) REAL*8 V(3,1) C V(1,1) = -DSIN(D)*DCOS(A) V(2,1) = -DSIN(D)*DSIN(A) V(3,1) = +DCOS(D) RETURN END C SUBROUTINE DLDA (A,D,V) C C BERECHNUNG DES NACH ALPHA ABGELEITETEN VEKTORS DER RICHTUNGSKOS. C C EINGABE : A, D = ALPHA, DELTA (IM BOGENMASS) C AUSGABE : V = VEKTOR DER ABGELETETEN RICHTUNGSKOSINUSSE C C IMPLICIT REAL*8 (A-H,O-Z) REAL*8 V(3,1) C V(1,1) = -DCOS(D)*DSIN(A) V(2,1) = +DCOS(D)*DCOS(A) V(3,1) = 0.D0 RETURN END C 1. Name: Subroutine BEJD U. Bastian, Okt. 1984 C C 2. Zweck: Besselsche Epochen in Julianische Daten umwandeln. C C 3. Aufruf: C CALL BEJD(BE,TJD) C C 4. Parameter: C C BE Input Besselsche Epoche (Jahre) Real*8 C C TJD Output Julianisches (Ephemeris) Datum (Tage) Real*8 C C SUBROUTINE BEJD(BE,TJD) IMPLICIT REAL*8 (A-Z) C TJD=(BE-1900.D0)*365.242198781D0+2415020.31352D0 RETURN END C******************************************************************** SUBROUTINE RADCHA(RADX,C2,CWERT,IER,LUNERR) IMPLICIT DOUBLE PRECISION (A-H,O-Z) CHARACTER*(*) C2,CWERT CHARACTER*10 CBT CHARACTER*1 CHR(20),CZ(4) CHARACTER*2 CH,CM,CS,CD,CV DATA CZ/' ','1','2','3'/ ILE1= LEN(C2) ILE2= LEN(CWERT) RAD=RADX IE1=0 INULF=0 DO 99 I=ILE1,1,-1 IF (C2(I:I).NE.' ') GOTO 98 99 CONTINUE IER=90 GOTO 900 98 LAE=I IF (C2(I:I).NE.'F') GOTO 97 LAE=I-1 INULF=1 97 CONTINUE IF (C2(ILE1:ILE1).EQ.'F') ILE1 = ILE1 - 1 IF(ILE1.NE.ILE2) THEN IER = 93 GOTO 900 END IF DO 96 I=1,ILE2 CWERT(I:I)=' ' CHR(I)=' ' 96 CONTINUE DO 95 I=1,LAE CHR(I)=C2(I:I) 95 CONTINUE RA=1. PI = 3.1415926535897932D0 IER=0 IVA=1 IH=0 IRUND=0 K1=0 IF (RAD.LT.0.D0) IVA=-1 RAD = RAD * IVA IF (CHR(2).EQ.'H') IH=1 IF (CHR(1).EQ.'T') IH=1 IF (CHR(2).EQ.'T') IH=1 IF (IH.EQ.1) GOTO 50 C******************************************************************** C C C C******************************************************************** R1=RAD*180.D0/PI RD=DINT(R1) RDB=R1-RD R2=RDB*60.D0 RM=DINT(R2) RMB=R2-RM R3=RMB*60.D0 RS=DINT(R3) RSB=R3-RS IF (RS.LT.60.D0) GOTO 41 RS=RS-60.D0 RM=RM+1.D0 41 IF (RM.LT.60.D0) GOTO 42 RM=RM-60.D0 RD=RD+1.D0 42 IF (CHR(2).EQ.'M') GOTO 60 IF (CHR(2).EQ.'A') GOTO 60 IF (CHR(2).EQ.'L' .AND. RD.LT.360.D0) GOTO 60 IF (CHR(2).EQ.'D' .AND. RD.LT.90.D0) GOTO 60 IF (CHR(2).EQ.'L') RD=RD-360.D0 IF (CHR(2).EQ.'D') RD=RD-90.D0 GOTO 42 C******************************************************************** C C C C******************************************************************** 50 R1=RAD*180.D0/(PI*15.D0) RH=DINT(R1) RHB=R1-RH R2=RHB*60.D0 RM=DINT(R2) RMB=R2-RM R3=RMB*60.D0 RS=DINT(R3) RSB=R3-RS IF (RS.LT.60.D0) GOTO 51 RS=RS-60.D0 RM=RM+1.D0 51 IF (RM.LT.60.D0) GOTO 52 RM=RM-60.D0 RH=RH+1.D0 52 IF (RH.LT.24.D0) GOTO 60 RH=RH-24.D0 GOTO 52 C******************************************************************** C C C C******************************************************************** 60 IF (CHR(1).EQ.'H') GOTO 100 IF (CHR(2).EQ.'D') GOTO 200 IF (CHR(2).EQ.'L') GOTO 500 IF (CHR(2).EQ.'T') GOTO 600 IF (CHR(2).EQ.'A') GOTO 700 IF (CHR(2).EQ.'M') GOTO 700 C******************************************************************** C C C C******************************************************************** 100 IE1=0 IF (INULF.EQ.1) GOTO 180 CALL INPUT1(RH,CH,IE) IE1=IE1+IE CALL INPUT1(RM,CM,IE) IE1=IE1+IE CALL INPUT1(RS,CS,IE) IE1=IE1+IE IF (IE1.EQ.0) GOTO 190 IER=101 GOTO 900 180 CALL INPUT2(RH,CH,IE) IE1=IE1+IE CALL INPUT2(RM,CM,IE) IE1=IE1+IE CALL INPUT2(RS,CS,IE) IE1=IE1+IE IF (IE1.EQ.0) GOTO 190 IER=102 GOTO 900 190 CWERT(1:2)=CH J=2 101 J=J+1 IF (J.GT.LAE) GOTO 401 IF (CHR(J).NE.' ') GOTO 102 CWERT(J:J)=' ' GOTO 101 102 IF (CHR(J).NE.'M') GOTO 110 C CWERT=CWERT//CM J1=J J2=J+1 CWERT(J1:J2)=CM J=J+1 103 J=J+1 IF (J.GT.LAE) GOTO 402 IF (CHR(J).NE.' ') GOTO 104 CWERT(J:J)=' ' GOTO 103 104 IF (CHR(J).NE.'S') GOTO 120 J1=J J2=J+1 CWERT(J1:J2)=CS J=J+1 C J=J+2 IF (J.GT.LAE) GOTO 403 105 J=J+1 IF (J.GT.LAE) GOTO 403 IF (CHR(J).NE.' ') GOTO 130 CWERT(J:J)=' ' GOTO 105 110 IF (CHR(J).NE.'.') GOTO 112 CWERT(J:J)='.' 111 J=J+1 IF (J.GT.LAE) GOTO 401 IF (CHR(J).NE.' ') GOTO 112 CWERT(J:J)=' ' GOTO 111 112 IF (CHR(J).NE.'H') GOTO 113 BT=RHB K1=1 GOTO 350 113 IF (CHR(J).NE.' ') GOTO 114 CWERT(J:J)=' ' J=J+1 IF (J.GT.LAE) GOTO 401 GOTO 112 114 IER=111 GOTO 900 120 IF (CHR(J).NE.'.') GOTO 122 CWERT(J:J)='.' 121 J=J+1 IF (J.GT.LAE) GOTO 402 IF (CHR(J).NE.' ') GOTO 122 CWERT(J:J)=' ' GOTO 121 122 IF (CHR(J).NE.'M') GOTO 123 BT=RMB K1=2 GOTO 350 123 IF (CHR(J).NE.' ') GOTO 124 CWERT(J:J)=' ' J=J+1 IF (J.GT.LAE) GOTO 402 GOTO 122 124 IER=121 GOTO 900 130 IF (CHR(J).NE.'.') GOTO 132 CWERT(J:J)='.' 131 J=J+1 IF (J.GT.LAE) GOTO 403 IF (CHR(J).NE.' ') GOTO 132 CWERT(J:J)=' ' GOTO 131 132 IF (CHR(J).NE.'S') GOTO 133 BT=RSB K1=3 GOTO 350 133 IF (CHR(J).NE.' ') GOTO 134 CWERT(J:J)=' ' J=J+1 IF (J.GT.LAE) GOTO 403 GOTO 132 134 IER=131 GOTO 900 C******************************************************************** C C C C******************************************************************** 200 IE1=0 IF (INULF.EQ.1) GOTO 280 CALL INPUT1(RD,CD,IE) IE1=IE1+IE CALL INPUT1(RM,CM,IE) IE1=IE1+IE CALL INPUT1(RS,CS,IE) IE1=IE1+IE IF (IE1.EQ.0) GOTO 290 IER=201 GOTO 900 280 CALL INPUT2(RD,CD,IE) IE1=IE1+IE CALL INPUT2(RM,CM,IE) IE1=IE1+IE CALL INPUT2(RS,CS,IE) IE1=IE1+IE IF (IE1.EQ.0) GOTO 290 IER=202 GOTO 900 290 CWERT(2:3)=CD J=3 201 J=J+1 IF (J.GT.LAE) GOTO 411 IF (CHR(J).NE.' ') GOTO 202 CWERT(J:J)=' ' GOTO 201 202 IF (CHR(J).NE.'M') GOTO 210 J1=J J2=J+1 CWERT(J1:J2)=CM J=J+1 203 J=J+1 IF (J.GT.LAE) GOTO 412 IF (CHR(J).NE.' ') GOTO 204 CWERT(J:J)=' ' GOTO 203 204 IF (CHR(J).NE.'S') GOTO 220 J1=J J2=J+1 CWERT(J1:J2)=CS J=J+1 IF (J.GT.LAE) GOTO 413 205 J=J+1 IF (J.GT.LAE) GOTO 413 IF (CHR(J).NE.' ') GOTO 230 CWERT(J:J)=' ' GOTO 205 210 IF (CHR(J).NE.'.') GOTO 212 CWERT(J:J)='.' 211 J=J+1 IF (J.GT.LAE) GOTO 411 IF (CHR(J).NE.' ') GOTO 212 CWERT(J:J)=' ' GOTO 211 212 IF (CHR(J).NE.'D') GOTO 213 BT=RDB K1=4 GOTO 350 213 IF (CHR(J).NE.' ') GOTO 214 CWERT(J:J)=' ' J=J+1 IF (J.GT.LAE) GOTO 411 GOTO 212 214 IER=211 GOTO 900 220 IF (CHR(J).NE.'.') GOTO 222 CWERT(J:J)='.' 221 J=J+1 IF (J.GT.LAE) GOTO 412 IF (CHR(J).NE.' ') GOTO 222 CWERT(J:J)=' ' GOTO 221 222 IF (CHR(J).NE.'M') GOTO 223 BT=RMB K1=5 GOTO 350 223 IF (CHR(J).NE.' ') GOTO 224 CWERT(J:J)=' ' J=J+1 IF (J.GT.LAE) GOTO 412 GOTO 222 224 IER=221 GOTO 900 230 IF (CHR(J).NE.'.') GOTO 232 CWERT(J:J)='.' 231 J=J+1 IF (J.GT.LAE) GOTO 413 IF (CHR(J).NE.' ') GOTO 232 CWERT(J:J)=' ' GOTO 231 232 IF (CHR(J).NE.'S') GOTO 233 BT=RSB K1=6 GOTO 350 233 IF (CHR(J).NE.' ') GOTO 234 CWERT(J:J)=' ' J=J+1 IF (J.GT.LAE) GOTO 413 GOTO 232 234 IER=231 GOTO 900 C******************************************************************** C C C C******************************************************************** 500 IE1=0 IF (RD.LT.360.D0) GOTO 591 RD=RD-360.D0 GOTO 500 591 RD1=DINT(RD/100.D0) RD2=RD-RD1*100.D0 IF (INULF.EQ.1) GOTO 580 CALL INPUT1(RD2,CD,IE) IE1=IE1+IE CALL INPUT1(RM,CM,IE) IE1=IE1+IE CALL INPUT1(RS,CS,IE) IE1=IE1+IE IF (IE1.EQ.0) GOTO 590 IER=301 GOTO 900 580 CALL INPUT2(RD2,CD,IE) IE1=IE1+IE CALL INPUT2(RM,CM,IE) IE1=IE1+IE CALL INPUT2(RS,CS,IE) IE1=IE1+IE IF (IE1.EQ.0) GOTO 590 IER=302 GOTO 900 590 CONTINUE DO 599 I=0,3 IF (RD1.EQ.DBLE(I)) GOTO 598 599 CONTINUE IER=301 GOTO 900 598 CWERT(2:2)=CZ(I+1) CWERT(3:4)=CD IF (CWERT(2:2).NE.' '.AND.CWERT(3:3).EQ.' ') CWERT(3:3)='0' IF (INULF.NE.1) GOTO 597 IF (CWERT(2:2).EQ.' ') CWERT(2:2)='0' IF (CWERT(3:3).EQ.' ') CWERT(3:3)='0' 597 J=4 501 J=J+1 IF (J.GT.LAE) GOTO 421 IF (CHR(J).NE.' ') GOTO 502 CWERT(J:J)=' ' GOTO 501 502 IF (CHR(J).NE.'M') GOTO 510 J1=J J2=J+1 CWERT(J1:J2)=CM J=J+1 503 J=J+1 IF (J.GT.LAE) GOTO 422 IF (CHR(J).NE.' ') GOTO 504 CWERT(J:J)=' ' GOTO 503 504 IF (CHR(J).NE.'S') GOTO 520 J1=J J2=J+1 CWERT(J1:J2)=CS J=J+1 IF (J.GT.LAE) GOTO 423 505 J=J+1 IF (J.GT.LAE) GOTO 423 IF (CHR(J).NE.' ') GOTO 530 CWERT(J:J)=' ' GOTO 505 510 IF (CHR(J).NE.'.') GOTO 512 CWERT(J:J)='.' 511 J=J+1 IF (J.GT.LAE) GOTO 421 IF (CHR(J).NE.' ') GOTO 512 CWERT(J:J)=' ' GOTO 511 512 IF (CHR(J).NE.'L') GOTO 513 BT=RDB K1=7 GOTO 350 513 IF (CHR(J).NE.' ') GOTO 514 CWERT(J:J)=' ' J=J+1 IF (J.GT.LAE) GOTO 421 GOTO 512 514 IER=311 GOTO 900 520 IF (CHR(J).NE.'.') GOTO 522 CWERT(J:J)='.' 521 J=J+1 IF (J.GT.LAE) GOTO 422 IF (CHR(J).NE.' ') GOTO 522 CWERT(J:J)=' ' GOTO 521 522 IF (CHR(J).NE.'M') GOTO 523 BT=RMB K1=8 GOTO 350 523 IF (CHR(J).NE.' ') GOTO 524 CWERT(J:J)=' ' J=J+1 IF (J.GT.LAE) GOTO 422 GOTO 522 524 IER=321 GOTO 900 530 IF (CHR(J).NE.'.') GOTO 532 CWERT(J:J)='.' 531 J=J+1 IF (J.GT.LAE) GOTO 423 IF (CHR(J).NE.' ') GOTO 532 CWERT(J:J)=' ' GOTO 531 532 IF (CHR(J).NE.'S') GOTO 533 BT=RSB K1=9 GOTO 350 533 IF (CHR(J).NE.' ') GOTO 534 CWERT(J:J)=' ' J=J+1 IF (J.GT.LAE) GOTO 423 GOTO 532 534 IER=331 GOTO 900 C******************************************************************** C C C C******************************************************************** 600 IE1=0 IF (INULF.EQ.1) GOTO 680 CALL INPUT1(RM,CM,IE) IE1=IE1+IE CALL INPUT1(RS,CS,IE) IE1=IE1+IE IF (IE1.EQ.0) GOTO 690 IER=101 GOTO 900 680 CALL INPUT2(RM,CM,IE) IE1=IE1+IE CALL INPUT2(RS,CS,IE) IE1=IE1+IE IF (IE1.EQ.0) GOTO 690 IER=102 GOTO 900 690 CWERT(2:3)=CM J=3 603 J=J+1 IF (J.GT.LAE) GOTO 432 IF (CHR(J).NE.' ') GOTO 604 CWERT(J:J)=' ' GOTO 603 604 IF (CHR(J).NE.'S') GOTO 620 J1=J J2=J+1 CWERT(J1:J2)=CS J=J+1 IF (J.GT.LAE) GOTO 433 605 J=J+1 IF (J.GT.LAE) GOTO 433 IF (CHR(J).NE.' ') GOTO 630 CWERT(J:J)=' ' GOTO 605 620 IF (CHR(J).NE.'.') GOTO 622 CWERT(J:J)='.' 621 J=J+1 IF (J.GT.LAE) GOTO 432 IF (CHR(J).NE.' ') GOTO 622 CWERT(J:J)=' ' GOTO 621 622 IF (CHR(J).NE.'M') GOTO 623 BT=RMB K1=11 GOTO 350 623 IF (CHR(J).NE.' ') GOTO 624 CWERT(J:J)=' ' J=J+1 IF (J.GT.LAE) GOTO 432 GOTO 622 624 IER=121 GOTO 900 630 IF (CHR(J).NE.'.') GOTO 632 CWERT(J:J)='.' 631 J=J+1 IF (J.GT.LAE) GOTO 433 IF (CHR(J).NE.' ') GOTO 632 CWERT(J:J)=' ' GOTO 631 632 IF (CHR(J).NE.'S') GOTO 633 BT=RSB K1=12 GOTO 350 633 IF (CHR(J).NE.' ') GOTO 634 CWERT(J:J)=' ' J=J+1 IF (J.GT.LAE) GOTO 433 GOTO 632 634 IER=131 GOTO 900 C******************************************************************** C C C C C******************************************************************** 700 IE1=0 IF (INULF.EQ.1) GOTO 780 CALL INPUT1(RM,CM,IE) IE1=IE1+IE CALL INPUT1(RS,CS,IE) IE1=IE1+IE IF (IE1.EQ.0) GOTO 790 IER=201 GOTO 900 780 CALL INPUT2(RM,CM,IE) IE1=IE1+IE CALL INPUT2(RS,CS,IE) IE1=IE1+IE IF (IE1.EQ.0) GOTO 790 IER=202 GOTO 900 790 CWERT(2:3)=CM J=3 703 J=J+1 IF (J.GT.LAE) GOTO 442 IF (CHR(J).NE.' ') GOTO 704 CWERT(J:J)=' ' GOTO 703 704 IF (CHR(J).NE.'S') GOTO 720 J1=J J2=J+1 CWERT(J1:J2)=CS J=J+1 IF (J.GT.LAE) GOTO 443 705 J=J+1 IF (J.GT.LAE) GOTO 443 IF (CHR(J).NE.' ') GOTO 730 CWERT(J:J)=' ' GOTO 705 720 IF (CHR(J).NE.'.') GOTO 722 CWERT(J:J)='.' 721 J=J+1 IF (J.GT.LAE) GOTO 442 IF (CHR(J).NE.' ') GOTO 722 CWERT(J:J)=' ' GOTO 721 722 IF (CHR(J).NE.'M') GOTO 723 BT=RMB K1=13 GOTO 350 723 IF (CHR(J).NE.' ') GOTO 724 CWERT(J:J)=' ' J=J+1 IF (J.GT.LAE) GOTO 442 GOTO 722 724 IER=221 GOTO 900 730 IF (CHR(J).NE.'.') GOTO 732 CWERT(J:J)='.' 731 J=J+1 IF (J.GT.LAE) GOTO 443 IF (CHR(J).NE.' ') GOTO 732 CWERT(J:J)=' ' GOTO 731 732 IF (CHR(J).NE.'S') GOTO 733 BT=RSB K1=15 GOTO 350 733 IF (CHR(J).NE.' ') GOTO 734 CWERT(J:J)=' ' J=J+1 IF (J.GT.LAE) GOTO 443 GOTO 732 734 IER=231 GOTO 900 C******************************************************************** C C C C******************************************************************** 350 IF (IRUND.NE.0) GOTO 380 K=LAE-J+1 A=DBLE(K) BT=DINT(BT*10.**A+0.5D0) AT=BT/(10.D0**A) IF (AT.LT.1.D0) GOTO 380 BT=BT-10.D0**(A+1.D0) IRUND=1 IF (K1.EQ.1) GOTO 361 IF (K1.EQ.4) GOTO 364 IF (K1.EQ.7) GOTO 367 IF (K1.EQ.2) GOTO 362 IF (K1.EQ.5) GOTO 365 IF (K1.EQ.8) GOTO 368 IF (K1.EQ.3) GOTO 363 IF (K1.EQ.6) GOTO 366 IF (K1.EQ.9) GOTO 369 IF (K1.EQ.11) GOTO 371 IF (K1.EQ.12) GOTO 372 IF (K1.EQ.14) GOTO 374 IF (K1.EQ.15) GOTO 375 361 RH=RH+1 IF (RH.EQ.24.) RH=0.D0 IF (RH.EQ.25.) RH=1.D0 GOTO 100 362 RM=RM+1 IF (RM.LT.60.D0) GOTO 100 RM=RM-60.D0 GOTO 361 363 RS=RS+1 IF (RS.LT.60.D0) GOTO 100 RS=RS-60.D0 GOTO 362 364 RD=RD+1 GOTO 200 365 RM=RM+1 IF (RM.LT.60.D0) GOTO 200 RM=RM-60.D0 GOTO 364 366 RS=RS+1 IF (RS.LT.60.D0) GOTO 200 RS=RS-60.D0 GOTO 365 367 RD=RD+1 IF (RD.GE.360.D0) RD=RD-360.D0 GOTO 500 368 RM=RM+1 IF (RM.LT.60.D0) GOTO 500 RM=RM-60.D0 GOTO 367 369 RS=RS+1 IF (RS.LT.60.D0) GOTO 500 RS=RS-60.D0 GOTO 368 371 RM=RM+1 IF (RM.LT.60.D0) GOTO 600 RM=RM-60.D0 GOTO 600 372 RS=RS+1 IF (RS.LT.60.D0) GOTO 600 RS=RS-60.D0 GOTO 371 374 RM=RM+1 IF (RM.LT.60.D0) GOTO 700 RM=RM-60.D0 GOTO 700 375 RS=RS+1 IF (RS.LT.60.D0) GOTO 700 RS=RS-60.D0 GOTO 374 C******************************************************************** C C C C******************************************************************** 380 CALL INPUT3(BT,K,CBT,IE) J1=J J2=LAE CWERT(J1:J2)=CBT GOTO 900 C******************************************************************** C C C C******************************************************************** 401 IF (IRUND.EQ.1) GOTO 900 IRUND=1 RHB1=RHB+0.5D0 IF (RHB1.LT.1.0D0) GOTO 900 GOTO 361 402 IF (IRUND.EQ.1) GOTO 900 IRUND=1 RMB1=RMB+0.5D0 IF (RMB1.LT.1.0D0) GOTO 900 GOTO 362 403 IF (IRUND.EQ.1) GOTO 900 IRUND=1 RSB1=RSB+0.5D0 IF (RSB1.LT.1.0D0) GOTO 900 GOTO 363 411 IF (IRUND.EQ.1) GOTO 900 IRUND=1 RDB1=RDB+0.5D0 IF (RDB1.LT.1.0D0) GOTO 900 GOTO 364 412 IF (IRUND.EQ.1) GOTO 900 IRUND=1 RMB1=RMB+0.5D0 IF (RMB1.LT.1.0D0) GOTO 900 GOTO 365 413 IF (IRUND.EQ.1) GOTO 900 IRUND=1 RSB1=RSB+0.5D0 IF (RSB1.LT.1.0D0) GOTO 900 GOTO 366 421 IF (IRUND.EQ.1) GOTO 900 IRUND=1 RDB1=RDB+0.5D0 IF (RDB1.LT.1.0D0) GOTO 900 GOTO 367 422 IF (IRUND.EQ.1) GOTO 900 IRUND=1 RMB1=RMB+0.5D0 IF (RMB1.LT.1.0D0) GOTO 900 GOTO 368 423 IF (IRUND.EQ.1) GOTO 900 IRUND=1 RSB1=RSB+0.5D0 IF (RSB1.LT.1.0D0) GOTO 900 GOTO 369 432 IF (IRUND.EQ.1) GOTO 900 IRUND=1 RMB1=RMB+0.5D0 IF (RMB1.LT.1.0D0) GOTO 900 GOTO 371 433 IF (IRUND.EQ.1) GOTO 900 IRUND=1 RSB1=RSB+0.5D0 IF (RSB1.LT.1.0D0) GOTO 900 GOTO 372 442 IF (IRUND.EQ.1) GOTO 900 IRUND=1 RMB1=RMB+0.5D0 IF (RMB1.LT.1.0D0) GOTO 900 GOTO 374 443 IF (IRUND.EQ.1) GOTO 900 IRUND=1 RSB1=RSB+0.5D0 IF (RSB1.LT.1.0D0) GOTO 900 GOTO 375 C******************************************************************** C C C C******************************************************************** 900 CONTINUE IF (IER.EQ.0) GOTO 910 WRITE(LUNERR,1111) RAD,C2,IER 1111 FORMAT (1X, 'FEHLER AUFGETRETEN BEI UMRECHNUNG VON'/ *1X,F11.8,' IN ',A/1X,'FEHLERMELDUNG = ',I3//) RAD = 0.D0 GOTO 920 910 CONTINUE CV='+' IF (CHR(1).EQ.'N') CV=' ' IF (IVA.EQ.-1) CV='-' IF (IH.NE.1) CWERT(1:1)=CV IF (IH.EQ.1 .AND. CHR(2).EQ.'T') CWERT(1:1)=CV 920 CONTINUE 999 CONTINUE RETURN END C******************************************************************** C C C C******************************************************************** SUBROUTINE INPUT1(Z,C,IE) IMPLICIT DOUBLE PRECISION (A-H,O-Z) CHARACTER*2 IZI(100),C DATA IZI/' 0',' 1',' 2',' 3',' 4',' 5',' 6',' 7',' 8',' 9', *'10','11','12','13','14','15','16','17','18','19', *'20','21','22','23','24','25','26','27','28','29', *'30','31','32','33','34','35','36','37','38','39', *'40','41','42','43','44','45','46','47','48','49', *'50','51','52','53','54','55','56','57','58','59', *'60','61','62','63','64','65','66','67','68','69', *'70','71','72','73','74','75','76','77','78','79', *'80','81','82','83','84','85','86','87','88','89', *'90','91','92','93','94','95','96','97','98','99'/ IE = 0 DO 10 I=0,99 IF (Z.EQ.DBLE(I)) C=IZI(I+1) 10 CONTINUE RETURN END C******************************************************************** C C C C******************************************************************** SUBROUTINE INPUT2(Z,C,IE) IMPLICIT DOUBLE PRECISION (A-H,O-Z) CHARACTER*2 IZI(100),C DATA IZI/'00','01','02','03','04','05','06','07','08','09', *'10','11','12','13','14','15','16','17','18','19', *'20','21','22','23','24','25','26','27','28','29', *'30','31','32','33','34','35','36','37','38','39', *'40','41','42','43','44','45','46','47','48','49', *'50','51','52','53','54','55','56','57','58','59', *'60','61','62','63','64','65','66','67','68','69', *'70','71','72','73','74','75','76','77','78','79', *'80','81','82','83','84','85','86','87','88','89', *'90','91','92','93','94','95','96','97','98','99'/ IE = 0 DO 10 I=0,99 IF (Z.EQ.DBLE(I)) C=IZI(I+1) 10 CONTINUE RETURN END C******************************************************************** C C C******************************************************************** SUBROUTINE INPUT3(B,L,CC,IE) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION Z(9) CHARACTER*10 CC CHARACTER*1 IZ(10),IZI(10) DATA IZ/'0','1','2','3','4','5','6','7','8','9'/ DO 99 I=1,10 CC(I:I)=' ' 99 CONTINUE IF (L.GT.9) GOTO 90 DO 10 I=1,9 Z(I)=0 IZI(I)=' ' 10 CONTINUE GOTO(1,2,3,4,5,6,7,8,9) L 9 Z(9)=DINT(B/10.**8.D0) 8 Z(8)=DINT((B-Z(9)*10.D0**8.D0)/10.**7.D0) 7 Z(7)=DINT((B-Z(9)*10.D0**8.D0-Z(8)*10.D0**7.D0)/10.**6.D0) 6 Z(6)=DINT((B-Z(9)*10.D0**8.D0-Z(8)*10.D0**7.D0- F Z(7)*10.D0**6.D0)/10.**5.D0) 5 Z(5)=DINT((B-Z(9)*10.D0**8.D0-Z(8)*10.D0**7.D0- F Z(7)*10.D0**6.D0-Z(6)*10.D0**5.D0)/10.**4.D0) 4 Z(4)=DINT((B-Z(9)*10.D0**8.D0-Z(8)*10.D0**7.D0- F Z(7)*10.D0**6.D0-Z(6)*10.D0**5.D0- F Z(5)*10.D0**4.D0)/1000.) 3 Z(3)=DINT((B-Z(9)*10.D0**8.D0-Z(8)*10.D0**7.D0- F Z(7)*10.D0**6.D0-Z(6)*10.D0**5.D0- F Z(5)*10.D0**4.D0-Z(4)*10.D0**3.D0)/100.) 2 Z(2)=DINT((B-Z(9)*10.D0**8.D0-Z(8)*10.D0**7.D0- F Z(7)*10.D0**6.D0-Z(6)*10.D0**5.D0- F Z(5)*10.D0**4.D0-Z(4)*10.D0**3.D0- F Z(3)*10.D0**2.D0)/10.) 1 Z(1)= DINT(B-Z(9)*10.D0**8.D0-Z(8)*10.D0**7.D0- F Z(7)*10.D0**6.D0-Z(6)*10.D0**5.D0- F Z(5)*10.D0**4.D0-Z(4)*10.D0**3.D0- F Z(3)*10.D0**2.D0-Z(2)*10.D0**1.D0) DO 12 I=0,9 DO 11 J=1,9 IF (Z(J).EQ.DBLE(I)) IZI(J)=IZ(I+1) 11 CONTINUE 12 CONTINUE CC(1:1)=IZI(L) DO 13 I=L-1,1,-1 L1=L-I+1 CC(L1:L1)=IZI(I) 13 CONTINUE GOTO 14 90 IE=1 14 CONTINUE RETURN C******************************************************************** C C C C******************************************************************** END C SUBROUTINE ANGLE(V,A,D) C C BERECHNUNG VON ALPHA, DELTA AUS DEM VEKTOR DER RICHTUNGSKOSIN. C C EINGABE : V = VEKTOR DER RICHTUNGSKOSINUSSE C AUSGABE : A,D = ALPHA, DELTA C C IMPLICIT REAL*8 (A-H,O-Z) REAL*8 V(3,1) C C PI = 3.141592653589793D0 C C C UEBERGANG ZUM EINHEITSVEKTOR. C DIES IST DANN NOTWENDIG, WENN VERSEHENTLICH NICHT DIE C RICHTUNGSKOSINUSSE SELBST, SONDERN EIN VEKTOR, DER EIN C VIELFACHES DAVON IST, EINGEGEBEN WURDE C CALL SCLPRD (V,V,X,3) X = DSQRT(X) DO 10 I = 1,3 10 V(I,1) = V(I,1)/X C D = ASIN(V(3,1)) A = DATAN2(V(2,1),V(1,1)) IF ( A.LT.0.D0 ) A = A + 2.D0*PI RETURN END C SUBROUTINE SCLPRD (A,B,R,M) C C SKALRAPRODUKT R = A.B ZWEIER VEKTOREN C C A = NAME DES 1. EINGABEVEKTORS (M-MAL-1 MATRIX) C B = NAME DES 2. EINGABEVEKTORS (M-MAL-1 MATRIX) C R = SKALAR-PRODUKT C M = ZAHL DER ZEILEN VON A = ZAHL DER ZEILEN VON B C C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(M,1),B(M,1) C R = 0.D0 C DO 10 I = 1,M 10 R = R + A(I,1)*B(I,1) RETURN END SUBROUTINE BECORA(L,M,MMP1,X,N,XMD,SD,RD,SXI,SXISQ,SXIXJ) C C SUBROUTINE ZUR SCHAETZUNG VON MITTELWERTEN UND STREUUNGEN, C SOWIE KORRELATIONSKOEFFIZIENTEN. C ERSETZT DIE WIRKLICH DUEMMLICHE IMSL-ROUTINE BECOR. C C PARAMETER: C C 1. INPUT-PARAMETER C C L INTEGER STEUERPARAMETER: C L=1: ES WIRD NOCH EINE DATENZEILE EINGEGEBEN. C DIE AUSGABEWERTE WERDEN NICHT BERECHNET. C L=2: ES WIRD KEINE DATENZEILE MEHR EINGEGEBEN. C DIE AUSGABEWERTE SOLLEN BERECHNET WERDEN. C L=SONST.: ES WIRD KEINE DATENZEILE MEHR EINGEGEBEN, C DIE ZAEHLER ETC. DER SUBROUTINE SOLLEN C INITIALISIERT WERDEN. C M INTEGER LAENGE DER VEKTOREN X,XMD,SD,SXI,SXISQ C ANZAHL DER VARIABLEN, FUER DIE DIE AUSGABEWERTE C BERECHNET WERDEN SOLLEN. C MMP1 INTEGER LAENGE DER VEKTOREN RD, SXIXJ C MINDESTENS M*(M+1)/2 C X REAL*8 VEKTOR DER MESSWERTE, FUER JEDE VARIABLE WIRD C BEI JEDEM AUFRUF DER SUBROUTINE GENAU EIN WERT C UEBERGEBEN. C WIRD NUR BENUTZT FUER L=1 C C 2. OUTPUT FUER L=1: C C N INTEGER ZAEHLER, DER DIE BEREITS EINGEGEBENEN C UND AKKUMULIERTEN DATENZEILEN ZAEHLT. C MIT ANDEREN WORTEN: NUMMER DER GERADE EINGEGEBENEN C DATENZEILE. C BEGINNT MIT 1 ZU LAUFEN BEIM ERSTEN AUFRUF C DER SUBROUTINE UND NACH JEDEM AUFRUF MIT L=2. C C N IST GLEICHZEITIG EIN FEHLERINDIKATOR: C FUER N=1 SIND ALLE VARIANZEN UND KOVARIANZEN AUF C NULL GESETZT. C C XMD REAL*8 VEKTOR DER LAENGE M C ENTHAELT DIE MITTELWERTE DER EINZELNEN VARIABLEN C C SD REAL*8 VEKTOR DER LAENGE M C ENTHAELT DIE STANDARDABWEICHUNGEN (ACHTUNG: NICHT C VARIANZEN) DER EINZELNEN VARIABLEN C C RD REAL*8 VEKTOR DER LAENGE MMP1 C ENTHAELT DIE KOVARIANZEN (ACHTUNG: NICHT C KORRELATIONSKOEFFIZIENTEN) DER VARIABLEN C IN SYMMETRIC STORAGE MODE (INDEX=I*(I-1)+J). C C SXI REAL*8 HILFSVEKTOR DER LAENGE M C ENTHAELT DIE SUMMEN UEBER X(I) C C SXISQ REAL*8 HILFSVEKTOR DER LAENGE M C ENTHAELT DIE SUMMEN UEBER X(I) MAL X(I) C C SXIXJ REAL*8 HILFSVEKTOR DER LAENGE MMP1 C ENTHAELT DIE SUMMEN UEBER X(I) MAL X(J) C IN SYMMETRIC STORAGE MODE C IMPLICIT REAL*8(A-H,O-Z) DIMENSION X(M),XMD(M),SD(M),RD(MMP1) DIMENSION SXI(M),SXISQ(M),SXIXJ(MMP1) C SAVE C C STATEMENT FUNCTION K(I,J)=I*(I-1)/2+J C DATA NTEST/0/ C C IF(NTEST.EQ.0.OR.L.GE.3) *THEN C C C INITIALISIEREN C N=0 DO 103 I=1,M SXI(I)=0.D0 SXISQ(I)=0.D0 SD(I)=0.D0 XMD(I)=0.D0 SD(I)=0.D0 DO 103 J=1,I SXIXJ(K(I,J))=0.D0 RD(K(I,J))=0.D0 103 CONTINUE C C END IF IF(L.EQ.1) *THEN C C C AKKUMULIEREN C N=N+1 DO 101 I=1,M SXI(I)=SXI(I)+X(I) SXISQ(I)=SXISQ(I)+X(I)*X(I) DO 101 J=1,I SXIXJ(K(I,J))=SXIXJ(K(I,J))+X(I)*X(J) 101 CONTINUE END IF C C AUSGABEWERTE BERECHNEN C IF(L.EQ.2) THEN DO 102 I=1,M XMD(I)=SXI(I)/N IF(N.GT.1) THEN SD(I)=DSQRT((SXISQ(I)-N*XMD(I)*XMD(I))/(N-1)) DO 104 J=1,I RD(K(I,J))=(SXIXJ(K(I,J))-N*XMD(I)*XMD(J))/(N-1) 104 CONTINUE END IF 102 CONTINUE C C END IF C C NTEST=1 RETURN END C***********************************************************************00000100 C 00000200 SUBROUTINE AUSR(SLDIS,V,H,A,G,Q,ABCRMS,HH,IL,KDER, 00000300 *IUU,N,MM,ILL,IAU,IAU0,IV,IER) 00000400 C 00000500 C***********************************************************************00000600 C 00000700 C1. Name: Subroutine A U S R H.Bernstein, 1988, mit 00000800 C======== Aenderungen von U.Bastian 00000900 C 00001000 C2. Zweck: Ausreissertest bei Ausgleichungsproblemen 00001100 C========= (Vergleich der studentisierten Residuen nach der 00001200 C Ausgleichung mit den Varianzen dieser Residuen) 00001300 C 00001400 C3. Aufruf: 00001500 C========== 00001600 C 00001700 C CALL AUSR(SLDIS,V,H,A,G,Q,ABCRMS,HH,IL,KDER, 00001800 C * IUU,N,MM,ILL,IAU,IAU0,IV,IER) 00001900 C 00002000 C4. Parameter: 00002100 C============= 00002200 C 00002300 CSLDIS INPUT Real*8 Signifikanzniveau fuer den Test. 00002400 C Moegliche Werte: 0.90, 0.95, 0.99 00002500 C Fuer alle anderen Werte wird IER=1 gesetzt und 00002600 C alle Ausgabedaten bleiben undefiniert. 00002700 CV INPUT Real*8(IL) Vektor der Residuen nach der Ausgleichung 00002800 CH INPUT Real*8(IL) Arbeitsvektor der Laenge IL 00002900 CA INPUT Real*8(IL,KDER) Koeffizientenmatrix (design matrix) 00003000 C des Ausgleichungsproblems. A(I,J) ist die 00003100 C Ableitung der i-ten Messung nach der j-ten 00003200 C Unbekannten. 00003300 CG INPUT Real*8(IUU) Kovarianzmatrix der Unbekannten (d.h. 00003400 C Inverse der Normalgleichungsmatrix) in 00003500 C symmetric storage mode. 00003600 CQ INPUT Real*8(ILL) Gewichtsmatrix der Messungen 00003700 C in symmetric storage mode 00003800 CABCRMS INPUT Real*8(IL) Mittlere Fehler der einzelnen Messungen 00003900 CHH INPUT Real*8(IL) Arbeitsvektor der Laenge IL 00004000 CIL INPUT Integer Maximale Anzahl der Messungen 00004100 CKDER INPUT Integer Maximale Anzahl der Unbekannten 00004200 CIUU INPUT Integer Maximale Anzahl der Elemente von G 00004300 C IUU muss groesser oder gleich N*(N+1)/2 sein. 00004400 CN INPUT Integer Tatsaechliche Anzahl der Unbekannten, die im 00004500 C Ausgleichungsproblem behandelt wurden. Damit werden 00004600 C aus der groesser dimensionierten A-Matrix nur die ersten 00004700 C N Spalten verwendet. Fuer die uebrigen Vektoren und 00004800 C Matrizen analog. 00004900 CMM INPUT Tatsaechliche Anzahl der Messungen, die im Ausgleichungs- 00005000 C problem behandelt wurden. Damit werden aus der groesser 00005100 C dimensionierten A-Matrix nur die ersten MM Zeilen 00005200 C verwendet. Fuer die uebrigen Matrizen und Vektoren analog.00005300 CILL INPUT Integer Dimension von Q. ILL=IL*(IL+1)/2 00005400 CIAU INPUT/OUTPUT INTEGER GESAMTZAHL DER BISHER ENTDECKTEN AUSREISSER00005500 C MUSS VOR DEM ERSTEN AUFRUF AUF NULL GESETZT WERDEN. 00005600 CIAU0 OUTPUT INTEGER ANZAHL DER AUSREISSER, DIE BEIM GEGENWAERTIGEN 00005700 C Aufruf der Subroutine entdeckt wurden. 00005800 CIV INPUT/OUTPUT Integer(IL) Vektor der die laufenden Nummern der 00005900 C bisher entdeckten Ausreisser angibt. 00006000 C IV(I)=0 falls die i-te Messung kein Ausreisser 00006100 C IV(I)=1 falls die i-te Messung Ausreisser ist. 00006200 C Muss vor dem ersten Aufruf auf null gesetzt werden. 00006300 CIER OUTPUT Integer Fehlercode. 00006400 C IER=0 alles o.k. 00006500 C IER=1 Fehler. Erklaerung s. oben bei SLDIS 00006600 C 00006700 C5. Notwendige Statements im rufenden Programm: 00006800 C============================================== 00006900 C 00007000 C Dimensionierungen. Nullsetzen der IAU und IV vor dem jeweils ersten 00007100 C Aufruf (fuer jede neue Ausgleichung|). 00007200 C 00007300 C6. Notwendige Steuerkarten: Keine 00007400 C=========================== 00007500 C 00007600 C7. Weitere benoetigte Subroutines: Keine 00007700 C================================== 00007800 C 00007900 C8. Genaue Funktionsbeschreibung: 00008000 C================================ 00008100 C 00008200 CDie Subroutine erkennt bei jedem Aufruf hoechstens einen 00008300 CAusreisser. Will man alle Ausreisser aus einem Ausgleichungsfall 00008400 Centfernen, so muss man die Subroutine so lange immer wieder aufrufen, 00008500 Cbis kein Ausreisser mehr gefunden wird (IAU0=0). 00008600 CZwischen den einzelnen Aufrufen sollte das Ausgleichungsproblem 00008700 Cerneut geloest werden, unter Auslassung der vorher bereits erkannten 00008800 CAusreisser. Dies ist aber nicht unbedingt erforderlich, wenn man 00008900 Cauf die volle Empfindlichkeit und statistische Konsistenz des 00009000 CVerfahrens verzichten will. 00009100 C 00009200 CNaeheres zur Arbeitsweise und zum statistischen Prinzip der 00009300 CSubroutine siehe in 00009400 C 00009500 CH.Wolf, Ausgleichsrechnung, Duemmler-Verlag, Bonn 1975, S. 275 oder 00009600 CFoerstner, Statistical Concepts for Quality Control 00009700 C9. Fundort fuer die Subroutine: 00009800 C=============================== 00009900 C 00010000 C FORTRAN-QUELLTEXT IN S.S19.ALL.SUBPROG.FORT 00010100 C LOAD-MODULE IN S.S19.ALL.FVLOAD 00010200 C 00010220 C 00010300 C 00010400 C 00010500 IMPLICIT REAL*8(A-H,O-Z) 00010600 C 00010700 DIMENSION V(IL),IV(IL),H(IL),G(IUU),A(IL,KDER),HH(IL),ABCRMS(IL) 00010800 DIMENSION Q(ILL) 00010900 COMMON/FILES/ IEDIT,NOUT,NSTOP,NINP,NSYS 00011000 PARAMETER (IC1ALP=90, IC3ALP=99, IC2ALP=95) 00011100 PARAMETER (C1RL0=6.1D0, C2RL0=7.9D0, C3RL0=11.6D0) 00011200 C 00011300 ISL = NINT (100.D0*SLDIS) 00011400 IAU0=0 00011500 C 00011600 00011700 IF(ISL.EQ.IC1ALP) THEN 00011800 RL0=C1RL0 00011900 ELSE IF(ISL.EQ.IC2ALP) THEN 00012000 RL0=C2RL0 00012100 ELSE IF(ISL.EQ.IC3ALP) THEN 00012200 RL0=C3RL0 00012300 ELSE 00012400 IER=1 00012500 RETURN 00012600 END IF 00012700 C 00012800 DO 100 I=1,MM 00012900 SS=0.0 00013000 DO 200 L=1,N 00013100 S=0.0 00013200 DO 300 J=1,N 00013300 IF(J.GT.L) JL=L+(J*J-J)/2 00013400 IF(J.LE.L) JL=J+(L*L-L)/2 00013500 S=S+A(I,J)*G(JL) 00013600 300 CONTINUE 00013700 HH(L)=S 00013800 200 CONTINUE 00013900 DO 400 J=1,N 00014000 SS=SS+HH(J)*A(I,J) 00014100 400 CONTINUE 00014200 H(I)=ABCRMS(I)*SQRT(ABS(RL0*(1.0-SS*Q(I+(I*I-I)/2)))) 00014300 100 CONTINUE 00014400 C 00014500 DO 500 I=1,MM 00014600 HH(I)=0.0 00014700 500 CONTINUE 00014800 C 00014900 DO 600 I=1,MM 00015000 IF (IV(I).EQ.1) GOTO 600 00015100 IF (ABS( V(I) ).GT.H(I)) HH(I)=ABS(ABS( V(I) )-H(I)) 00015200 600 CONTINUE 00015300 C 00015400 K=1 00015500 DO 700 I=1,MM 00015600 IF (HH(I).GT.HH(K)) K=I 00015700 700 CONTINUE 00015800 C 00015900 IF (HH(K).NE.0.0) THEN 00016000 IAU=IAU+1 00016100 IAU0=IAU0+1 00016200 IV(K)=1 00016300 END IF 00016400 C 00016500 RETURN 00016600 END 00016700 SUBROUTINE LINF2(N,T,X,W,IV,IAU,MT,MX,M,X0,SMX,SM,S0) 00000100 C-----------------------------------------------------------------------00000200 C ES WIRD EINE AUSGLEICHSGERADE DURCH DIE PUNKTE 00000300 C T(I),X(I) MIT I = 1, ..., N 00000400 C BESTIMMT. 00000500 C 00000600 C 00000700 C EINGANGSPARAMETER: 00000800 C 00000900 C N = ZAHL DER PUNKTE 00001000 C T(I) = EPOCHE 00001100 C X(I) = BEOBACHTETER 'ORT' ZUR EPOCHE T(I) 00001200 C W(I) = GEWICHT DES ORTS X(I) 00001300 C 00001400 C AUSGANGSPARAMETER: 00001500 C 00001600 C MT = MITTLERE EPOCHE 00001700 C MX = MITTLERER ORT 00001800 C M = EIGENBEWEGUNG 00001900 C X0 = BERECHNETER ORT ZUR EPOCHE T = 0. 00002000 C SMX = MITTLERER FEHLER DES MITTLEREN ORTS 00002100 C SM = MITTLERER FEHLER DER EIGENBEWEGUNG 00002200 C-----------------------------------------------------------------------00002500 REAL*8 T(N),X(N),MT,MX,M,X0,SMX,SM,D,V,S0Q,S0,TAU, 00002600 @ SW,W(N),MXQ 00002700 DIMENSION IV(N) 00002710 MT=0.D0 00002800 SW=0.D0 00002900 MX=0.D0 00003000 MXQ=0.D0 00003100 DO 10 I=1,N 00003200 IF(IV(I).EQ.1) GOTO 10 00003210 MT=MT+T(I)*W(I) 00003300 SW=SW+W(I) 00003400 MX=MX+X(I)*W(I) 00003500 MXQ=MXQ+X(I)*X(I)*W(I) 00003600 10 CONTINUE 00003610 MT=MT/SW 00003700 MX=MX/SW 00003800 C-----------------------------------------------------------------------00003900 C MXQ = MITTLERES X(I)**2 00004000 C-----------------------------------------------------------------------00004100 MXQ=MXQ/SW 00004200 M=0.D0 00004300 D=0.D0 00004400 DO 20 I=1,N 00004500 IF(IV(I).EQ.1) GOTO 20 00004520 TAU=T(I)-MT 00004600 M=M+TAU*X(I)*W(I) 00004700 D=D+TAU*TAU*W(I) 00004800 20 CONTINUE 00004810 M=M/D 00004900 S0Q=0.D0 00005000 DO 30 I=1,N 00005100 IF(IV(I).EQ.1) GOTO 30 00005120 TAU=T(I)-MT 00005200 V=X(I)-MX-M*TAU 00005300 S0Q=S0Q+V*V*W(I) 00005400 30 CONTINUE 00005410 IF((N-IAU).LE.2) THEN 00005412 S0=0.D0 00005414 ELSE 00005416 S0Q=S0Q/(N-2-IAU) 00005500 S0=DSQRT(S0Q) 00005510 END IF 00005520 C-----------------------------------------------------------------------00005600 C S0 IST DER MITTLERE FEHLER DER GEWICHTSEINHEIT, FALLS N-IAU.GT.2 00005700 C-----------------------------------------------------------------------00005800 SMX=S0/DSQRT(SW) 00006000 SM=S0/DSQRT(D) 00006100 X0=MX-M*MT 00006200 RETURN 00006600 END 00006700 *DECK MDCHI C IMSL ROUTINE NAME - MDCHI C C----------------------------------------------------------------------- C C COMPUTER - CDC/SINGLE C C LATEST REVISION - JANUARY 1, 1978 C C PURPOSE - INVERSE CHI-SQUARED PROBABILITY DISTRIBUTION C FUNCTION C C USAGE - CALL MDCHI (P,DF,X,IER) C C ARGUMENTS P - INPUT PROBABILITY IN THE EXCLUSIVE RANGE C (0,1) C DF - INPUT NUMBER OF DEGREES OF FREEDOM. DF MUST BE C IN THE EXCLUSIVE RANGE (.5,200000.). C X - OUTPUT CHI-SQUARED VALUE, SUCH THAT A RANDOM C VARIABLE, DISTRIBUTED AS CHI-SQUARED WITH C DF DEGREES OF FREEDOM, WILL BE LESS THAN OR C EQUAL TO X WITH PROBABILITY P. C IER - ERROR PARAMETER. (OUTPUT) C TERMINAL ERROR C IER = 129 INDICATES THAT THE BOUNDS WHICH C ENCLOSED P COULD NOT BE FOUND WITHIN 20 C (NCT) ITERATIONS C IER = 130 INDICATES THAT AN ERROR OCCURRED C IN IMSL ROUTINE MDNRIS (P IS NOT A VALID C VALUE) C IER = 131 INDICATES THAT AN ERROR OCCURRED C IN IMSL ROUTINE MDCH C IER = 132 INDICATES THAT THE VALUE X COULD C NOT BE FOUND WITHIN 50 (ITMAX) ITERATIONS, C SO THAT THE ABSOLUTE VALUE OF P1-P WAS C LESS THAN OR EQUAL TO EPS. (P1 IS THE C CALCULATED PROBABILITY AT X, EPS = .0001) C C PRECISION/HARDWARE - SINGLE/ALL C C REQD. IMSL ROUTINES - H32/MDCH,MDNOR,MDNRIS,MERFI,MERRC=ERFC, C MGAMAD=DGAMMA,UERTST,UGETIO C - H36,H48,H60/MDCH,MDNOR,MDNRIS,MERFI, C MERRC=ERFC,MGAMA=GAMMA,UERTST,UGETIO C C NOTATION - INFORMATION ON SPECIAL NOTATION AND C CONVENTIONS IS AVAILABLE IN THE MANUAL C INTRODUCTION OR THROUGH IMSL ROUTINE UHELP C C COPYRIGHT - 1978 BY IMSL, INC. ALL RIGHTS RESERVED. C C WARRANTY - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN C APPLIED TO THIS CODE. NO OTHER WARRANTY, C EXPRESSED OR IMPLIED, IS APPLICABLE. C C----------------------------------------------------------------------- C SUBROUTINE MDCHI (P,DF,X,IER) C DATA EPS/.0001/,ITMAX/50/,NCT/20/,NSIG/5/ C FIRST EXECUTABLE STATEMENT IC = 0 IER = 0 C ESTIMATE STARTING X CALL MDNRIS (P,XP,IER) IF (IER .NE. 0) GO TO 80 DFF = .22222222222222 DFF = DFF/DF D = SQRT(DFF) X = DF*(1. -DFF + XP*D)**3 C IS THE CASE ASYMPTOTICALLY NORMAL IF (DF .GE. 40.) GO TO 9005 C FIND BOUNDS (IN X) WHICH ENCLOSE P NCNT = 0 IST = 0 ISW = 0 DX = X * .125 5 IF (X) 10,15,20 10 X = 0.0 DX = -DX GO TO 20 15 DX = .1 20 CALL MDCH (X,DF,P1,IER) DX = DX + DX IF (IER .NE. 0) GO TO 85 CSS = X NCNT = NCNT + 1 IF (NCNT .GT. NCT) GO TO 90 IF (P1 - P) 25,9005,30 25 X = X + DX ISW = 1 IF (IST .EQ. 0) GO TO 5 GO TO 35 30 X = X - DX IST = 1 IF (ISW .EQ. 0) GO TO 5 XR = CSS XL = X GO TO 40 35 XL = CSS XR = X C PREPARE FOR ITERATION TO FIND X 40 EPSP = 10.**(-NSIG) IF (XL .LT. 0.) XL = 0.0 CALL MDCH (XL,DF,P1,IER) IF (IER .NE. 0) GO TO 85 FXL = P1 - P CALL MDCH (XR,DF,P1,IER) IF (IER .NE. 0) GO TO 85 FXR = P1-P IF (FXL*FXR .NE. 0.0) GO TO 45 X = XR IF (FXL .EQ. 0.0) X = XL GO TO 9005 45 IF (DF .LE. 2. .OR. P .GT. .98) GO TO 50 C REGULA FALSI METHOD X = XL + FXL*(XR-XL)/(FXL-FXR) GO TO 55 C BISECTION METHOD 50 X = (XL+XR) * .5 55 CALL MDCH (X,DF,P1,IER) IF (IER .NE. 0) GO TO 85 FCS = P1-P IF (ABS(FCS) .GT. EPS) GO TO 60 GO TO 9005 60 IF (FCS * FXL .GT. 0.0) GO TO 65 XR = X FXR = FCS GO TO 70 65 XL = X FXL = FCS 70 IF (XR-XL .GT. EPSP*ABS(XR)) GO TO 75 GO TO 9005 75 IC = IC+1 IF (IC .LE. ITMAX) GO TO 45 IER = 132 GO TO 9000 C ERROR RETURNED FROM MDNRIS 80 IER = 130 GO TO 9000 C ERROR RETURNED FROM MDCH 85 IER = 131 GO TO 9000 90 IER = 129 9000 CONTINUE CALL UERTST (IER,6HMDCHI ) 9005 RETURN END *DECK MDNOR C IMSL ROUTINE NAME - MDNOR C C----------------------------------------------------------------------- C C COMPUTER - CDC/SINGLE C C LATEST REVISION - JANUARY 1, 1978 C C PURPOSE - NORMAL OR GAUSSIAN PROBABILITY DISTRIBUTION C FUNCTION C C USAGE - CALL MDNOR (Y,P) C C ARGUMENTS Y - INPUT VALUE AT WHICH FUNCTION IS TO BE C EVALUATED. C P - OUTPUT PROBABILITY THAT A RANDOM VARIABLE C HAVING A NORMAL (0,1) DISTRIBUTION WILL BE C LESS THAN OR EQUAL TO Y. C C PRECISION/HARDWARE - SINGLE/ALL C C REQD. IMSL ROUTINES - MERRC=ERFC C C NOTATION - INFORMATION ON SPECIAL NOTATION AND C CONVENTIONS IS AVAILABLE IN THE MANUAL C INTRODUCTION OR THROUGH IMSL ROUTINE UHELP C C COPYRIGHT - 1978 BY IMSL, INC. ALL RIGHTS RESERVED. C C WARRANTY - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN C APPLIED TO THIS CODE. NO OTHER WARRANTY, C EXPRESSED OR IMPLIED, IS APPLICABLE. C C----------------------------------------------------------------------- C SUBROUTINE MDNOR (Y,P) C SPECIFICATIONS FOR ARGUMENTS REAL P,Y C SPECIFICATIONS FOR LOCAL VARIABLES REAL SQR1D2 DATA SQR1D2/.70710678118655/ C FIRST EXECUTABLE STATEMENT P = .5 * ERFC(-Y*SQR1D2) RETURN END *DECK MERRC C IMSL ROUTINE NAME - MERRC=ERFC C C----------------------------------------------------------------------- C C COMPUTER - CDC/SINGLE C C LATEST REVISION - JUNE 1, 1981 C C PURPOSE - EVALUATE THE COMPLEMENTED ERROR FUNCTION C C USAGE - RESULT = ERFC(Y) C C ARGUMENTS Y - INPUT ARGUMENT OF THE COMPLEMENTED ERROR C FUNCTION. C ERFC - OUTPUT VALUE OF THE COMPLEMENTED ERROR C FUNCTION. C C PRECISION/HARDWARE - SINGLE/ALL C NOTE - ERFC MAY NOT BE SUPPLIED BY IMSL IF IT C RESIDES IN THE MATHEMATICAL SUBPROGRAM C LIBRARY SUPPLIED BY THE MANUFACTURER. C C REQD. IMSL ROUTINES - NONE REQUIRED C C NOTATION - INFORMATION ON SPECIAL NOTATION AND C CONVENTIONS IS AVAILABLE IN THE MANUAL C INTRODUCTION OR THROUGH IMSL ROUTINE UHELP C C COPYRIGHT - 1978 BY IMSL, INC. ALL RIGHTS RESERVED. C C WARRANTY - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN C APPLIED TO THIS CODE. NO OTHER WARRANTY, C EXPRESSED OR IMPLIED, IS APPLICABLE. C C----------------------------------------------------------------------- C REAL FUNCTION ERFC(Y) C SPECIFICATIONS FOR ARGUMENTS REAL Y C SPECIFICATIONS FOR LOCAL VARIABLES INTEGER ISW,I DIMENSION P(5),Q(3),P1(8),Q1(7),P2(5),Q2(4) REAL P,Q,P1,Q1,P2,Q2,XMIN,XLARGE,SSQPI,X, * RES,XSQ,XNUM,XDEN,XI C COEFFICIENTS FOR 0.0 .LE. Y .LT. C .477 DATA P(1)/-.44422647396874/, 1 P(2)/10.731707253648/, 2 P(3)/15.915606197771/, 3 P(4)/374.81624081284/, 4 P(5)/2.5612422994823E-02/ DATA Q(1)/17.903143558843/, 1 Q(2)/124.82892031581/, 2 Q(3)/332.17224470532/ C COEFFICIENTS FOR .477 .LE. Y C .LE. 4.0 DATA P1(1)/7.2117582508831/, 1 P1(2)/43.162227222057/, 2 P1(3)/152.98928504694/, 3 P1(4)/339.32081673434/, 4 P1(5)/451.91895371187/, 5 P1(6)/300.45926102016/, 6 P1(7)/-1.3686485738272E-07/, 7 P1(8)/.56419551747897/ DATA Q1(1)/77.000152935229/, 1 Q1(2)/277.58544474399/, 2 Q1(3)/638.98026446563/, 3 Q1(4)/931.35409485061/, 4 Q1(5)/790.95092532790/, 5 Q1(6)/300.45926095698/, 6 Q1(7)/12.782727319629/ C COEFFICIENTS FOR 4.0 .LT. Y DATA P2(1)/-.22695659353969/, 1 P2(2)/-4.9473091062325E-02/, 2 P2(3)/-2.9961070770354E-03/, 3 P2(4)/-2.2319245973418E-02/, 4 P2(5)/-2.7866130860965E-01/ DATA Q2(1)/1.0516751070679/, 1 Q2(2)/.19130892610783/, 2 Q2(3)/1.0620923052847E-02/, 3 Q2(4)/1.9873320181714/ C CONSTANTS DATA XMIN/1.0E-8/,XLARGE/5.6875E0/ C ERFC(XBIG) .APPROX. SETAP DATA XBIG/25.90625/ DATA SSQPI/.56418958354776/ C FIRST EXECUTABLE STATEMENT X = Y ISW = 1 IF (X.GE.0.0E0) GO TO 5 ISW = -1 X = -X 5 IF (X.LT..477E0) GO TO 10 IF (X.LE.4.0E0) GO TO 30 IF (ISW .GT. 0) GO TO 40 IF (X.LT.XLARGE) GO TO 45 RES = 2.0E0 GO TO 65 C ABS(Y) .LT. .477, EVALUATE C APPROXIMATION FOR ERFC 10 IF (X.LT.XMIN) GO TO 20 XSQ = X*X XNUM = P(5) DO 15 I = 1,4 XNUM = XNUM*XSQ+P(I) 15 CONTINUE XDEN = ((Q(1)+XSQ)*XSQ+Q(2))*XSQ+Q(3) RES = X*XNUM/XDEN GO TO 25 20 RES = X*P(4)/Q(3) 25 IF (ISW.EQ.-1) RES = -RES RES = 1.0E0-RES GO TO 65 C .477 .LE. ABS(Y) .LE. 4.0 C EVALUATE APPROXIMATION FOR ERFC 30 XSQ = X*X XNUM = P1(7)*X+P1(8) XDEN = X+Q1(7) DO 35 I=1,6 XNUM = XNUM*X+P1(I) XDEN = XDEN*X+Q1(I) 35 CONTINUE RES = XNUM/XDEN GO TO 55 C 4.0 .LT. ABS(Y), EVALUATE C MINIMAX APPROXIMATION FOR ERFC 40 IF (X.GT.XBIG) GO TO 60 45 XSQ = X*X XI = 1.0E0/XSQ XNUM = P2(4)*XI+P2(5) XDEN = XI+Q2(4) DO 50 I = 1,3 XNUM = XNUM*XI+P2(I) XDEN = XDEN*XI+Q2(I) 50 CONTINUE RES = (SSQPI+XI*XNUM/XDEN)/X 55 RES = RES*EXP(-XSQ) IF (ISW.EQ.-1) RES = 2.0E0-RES GO TO 65 60 RES = 0.0E0 65 ERFC = RES RETURN END *DECK MGAMA C IMSL ROUTINE NAME - MGAMA=GAMMA C C----------------------------------------------------------------------- C C COMPUTER - CDC/SINGLE C C LATEST REVISION - JUNE 1, 1981 C C PURPOSE - EVALUATE THE GAMMA FUNCTION C C USAGE - RESULT = GAMMA(X) C C ARGUMENTS X - INPUT ARGUMENT. C GAMMA IS SET TO MACHINE INFINITY, WITH THE C PROPER SIGN, WHENEVER C X IS ZERO, C X IS A NEGATIVE INTEGER, C ABS(X) .LE. XMIN, OR C ABS(X) .GE. XMAX. C XMIN IS OF THE ORDER OF 10**(-39) AND C XMAX IS AT LEAST 34.8. THE EXACT VALUES C OF XMIN AND XMAX MAY ALLOW LARGER RANGES C FOR X ON SOME COMPUTERS. C SEE THE PROGRAMMING NOTES IN THE MANUAL C FOR THE EXACT VALUES. C GAMMA - OUTPUT SINGLE PRECISION VALUE OF THE GAMMA C FUNCTION. C C PRECISION/HARDWARE - SINGLE/ALL C NOTE - GAMMA MAY NOT BE SUPPLIED BY IMSL IF C IT RESIDES IN THE MATHEMATICAL SUBPROGRAM C LIBRARY SUPPLIED BY THE MANUFACTURER. C C REQD. IMSL ROUTINES - UERTST,UGETIO C C NOTATION - INFORMATION ON SPECIAL NOTATION AND C CONVENTIONS IS AVAILABLE IN THE MANUAL C INTRODUCTION OR THROUGH IMSL ROUTINE UHELP C C REMARKS AN ERROR MESSAGE PRINTED BY UERTST FROM GAMMA SHOULD C BE INTERPRETED AS FOLLOWS C IER - ERROR INDICATOR C TERMINAL ERROR C IER = 129 INDICATES THAT THE ABSOLUTE VALUE C OF THE INPUT ARGUMENT X WAS SPECIFIED C GREATER THAN OR EQUAL TO XMAX. GAMMA C IS SET TO MACHINE INFINITY. C IER = 130 INDICATES THAT THE INPUT ARGUMENT C X WAS SPECIFIED AS ZERO OR A NEGATIVE C INTEGER OR THAT THE ABSOLUTE VALUE OF THE C INPUT ARGUMENT WAS LESS THAN OR EQUAL TO C XMIN. GAMMA IS SET TO MACHINE INFINITY C WITH THE PROPER SIGN FOR THE GAMMA C FUNCTION. IF X IS ZERO OR AN EVEN C NEGATIVE INTEGER, GAMMA HAS A NEGATIVE C SIGN. OTHERWISE IT HAS A POSITIVE SIGN. C C COPYRIGHT - 1978 BY IMSL, INC. ALL RIGHTS RESERVED. C C WARRANTY - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN C APPLIED TO THIS CODE. NO OTHER WARRANTY, C EXPRESSED OR IMPLIED, IS APPLICABLE. C C----------------------------------------------------------------------- C REAL FUNCTION GAMMA (X) C SPECIFICATIONS FOR ARGUMENTS REAL X C SPECIFICATIONS FOR LOCAL VARIABLES REAL P,Q,P4,BIG1,PI,XMIN,SIGN,Y,T,R,A,TOP, * DEN,B INTEGER I,IEND,IEND1,IEND2,IER,J LOGICAL MFLAG DIMENSION P(7),Q(6),P4(5) C COEFFICIENTS FOR MINIMAX C APPROXIMATION TO GAMMA(X), C 2.0 .LE. X .LE. 3.0 DATA P(1)/3.4109112397125E+01/, 1 P(2)/ -4.8341273405598E+01/, 2 P(3)/4.3005887829594E+02/, 3 P(4)/-5.5688734338586E+01/, 4 P(5)/2.0585220673644E+03/, 5 P(6)/7.7192407739800E-01/, 6 P(7)/-3.1721064346240E+00/ DATA Q(1)/2.4455138217658E+02/, 1 Q(2)/-1.0174768492818E+03/, 2 Q(3)/1.1615998272754E+03/, 3 Q(4)/2.0512896777972E+03/, 4 Q(5)/6.8080353498091E-01/, 5 Q(6)/-2.5386729086746E+01/ C COEFFICIENTS FOR MINIMAX C APPROXIMATION TO LN(GAMMA(X)), C 12.0 .LE. X DATA P4(1)/9.1893853320467E-01/, 1 P4(2)/8.3333333333267E-02/, 2 P4(3)/-2.7777776505151E-03/, 3 P4(4)/7.9358449768E-04/, 4 P4(5)/-5.82399983E-04/ DATA IEND/7/,IEND1/6/,IEND2/5/ DATA PI/3.1415926535898/ C GAMMA(XMIN) .APPROX. R1MACH(1) C GAMMA(BIG1) .APPROX. R1MACH(1) DATA XMIN/0.0/ DATA BIG1/177.803/ C FIRST EXECUTABLE STATEMENT IER = 0 MFLAG = .FALSE. T = X IF (ABS(T).GT.XMIN) GO TO 5 IER = 130 GAMMA = R1MACH(1) IF (T.LE.0.0) GAMMA = -R1MACH(1) GO TO 9000 5 IF (ABS(T).LT.BIG1) GO TO 10 IER = 129 GAMMA = R1MACH(1) GO TO 9000 10 IF (T.GT.0.0) GO TO 25 C ARGUMENT IS NEGATIVE MFLAG = .TRUE. T = -T R = AINT(T) SIGN = 1.0 IF (AMOD(R,2.0) .EQ. 0.0) SIGN = -1. R = T-R IF (R.NE.0.0) GO TO 20 IER = 130 GAMMA = R1MACH(1) IF (SIGN.EQ.-1.0) GAMMA = -R1MACH(1) GO TO 9000 C ARGUMENT IS NOT A NEGATIVE INTEGER 20 R = PI/SIN(R*PI)*SIGN T = T+1.0 C EVALUATE APPROXIMATION FOR GAMMA(T) C T .GT. XMIN 25 IF (T.GT.12.0) GO TO 60 I = T A = 1.0 IF (I.GT.2) GO TO 40 I = I+1 GO TO (30,35,50),I C 0.0 .LT. T .LT. 1.0 30 A = A/(T*(T+1.0)) T = T+2.0 GO TO 50 C 1.0 .LE. T .LT. 2.0 35 A = A/T T = T+1.0 GO TO 50 C 3.0 .LE. T .LE. 12.0 40 DO 45 J=3,I T = T-1.0 A = A*T 45 CONTINUE C 2.0 .LE. T .LE. 3.0 50 TOP = P(IEND1)*T+P(IEND) DEN = T+Q(IEND1) DO 55 J=1,IEND2 TOP = TOP*T+P(J) DEN = DEN*T+Q(J) 55 CONTINUE Y = (TOP/DEN)*A IF (MFLAG) Y = R/Y GAMMA = Y GO TO 9005 C T .GT. 12.0 60 TOP = ALOG(T) TOP = (T-1.5)*(TOP-1.0)+TOP-1.5 T = 1.0/T B = T*T Y = (((P4(5)*B+P4(4))*B+P4(3))*B+P4(2))*T+P4(1)+TOP Y = EXP(Y) IF (MFLAG) Y = R/Y GAMMA = Y GO TO 9005 9000 CONTINUE CALL UERTST(-IER,6H MGAMA) CALL UERTST(IER,6HGAMMA ) 9005 RETURN END *DECK UERTST C IMSL ROUTINE NAME - UERTST C C----------------------------------------------------------------------- C C COMPUTER - CDC/SINGLE C C LATEST REVISION - JANUARY 1, 1978 C C PURPOSE - PRINT A MESSAGE REFLECTING AN ERROR CONDITION C C USAGE - CALL UERTST (IER,NAME) C C ARGUMENTS IER - ERROR PARAMETER. (INPUT) C IER = I+J WHERE C I = 128 IMPLIES TERMINAL ERROR, C I = 64 IMPLIES WARNING WITH FIX, AND C I = 32 IMPLIES WARNING. C J = ERROR CODE RELEVANT TO CALLING C ROUTINE. C NAME - A SIX CHARACTER LITERAL STRING GIVING THE C NAME OF THE CALLING ROUTINE. (INPUT) C C PRECISION/HARDWARE - SINGLE/ALL C C REQD. IMSL ROUTINES - UGETIO C C NOTATION - INFORMATION ON SPECIAL NOTATION AND C CONVENTIONS IS AVAILABLE IN THE MANUAL C INTRODUCTION OR THROUGH IMSL ROUTINE UHELP C C REMARKS THE ERROR MESSAGE PRODUCED BY UERTST IS WRITTEN C ONTO THE STANDARD OUTPUT UNIT. THE OUTPUT UNIT C NUMBER CAN BE DETERMINED BY CALLING UGETIO AS C FOLLOWS.. CALL UGETIO(1,NIN,NOUT). C THE OUTPUT UNIT NUMBER CAN BE CHANGED BY CALLING C UGETIO AS FOLLOWS.. C NIN = 0 C NOUT = NEW OUTPUT UNIT NUMBER C CALL UGETIO(3,NIN,NOUT) C SEE THE UGETIO DOCUMENT FOR MORE DETAILS. C C COPYRIGHT - 1978 BY IMSL, INC. ALL RIGHTS RESERVED. C C WARRANTY - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN C APPLIED TO THIS CODE. NO OTHER WARRANTY, C EXPRESSED OR IMPLIED, IS APPLICABLE. C C----------------------------------------------------------------------- C SUBROUTINE UERTST (IER,NAME) C SPECIFICATIONS FOR ARGUMENTS INTEGER IER CHARACTER NAME*6 C SPECIFICATIONS FOR LOCAL VARIABLES CHARACTER*6 NAMSET,NAMEQ,IEQ*1 DATA NAMSET/'UERSET'/ DATA NAMEQ/' '/ C FIRST EXECUTABLE STATEMENT DATA LEVEL/4/,IEQDF/0/,IEQ/'='/ IF (IER.GT.999) GO TO 25 IF (IER.LT.-32) GO TO 55 IF (IER.LE.128) GO TO 5 IF (LEVEL.LT.1) GO TO 30 C PRINT TERMINAL MESSAGE CALL UGETIO(1,NIN,IOUNIT) IF (IEQDF.EQ.1) WRITE(IOUNIT,35) IER,NAMEQ,IEQ,NAME IF (IEQDF.EQ.0) WRITE(IOUNIT,35) IER,NAME GO TO 30 5 IF (IER.LE.64) GO TO 10 IF (LEVEL.LT.2) GO TO 30 C PRINT WARNING WITH FIX MESSAGE CALL UGETIO(1,NIN,IOUNIT) IF (IEQDF.EQ.1) WRITE(IOUNIT,40) IER,NAMEQ,IEQ,NAME IF (IEQDF.EQ.0) WRITE(IOUNIT,40) IER,NAME GO TO 30 10 IF (IER.LE.32) GO TO 15 C PRINT WARNING MESSAGE IF (LEVEL.LT.3) GO TO 30 CALL UGETIO(1,NIN,IOUNIT) IF (IEQDF.EQ.1) WRITE(IOUNIT,45) IER,NAMEQ,IEQ,NAME IF (IEQDF.EQ.0) WRITE(IOUNIT,45) IER,NAME GO TO 30 15 CONTINUE C CHECK FOR UERSET CALL IF (NAME.NE.NAMSET) GO TO 25 LEVOLD = LEVEL LEVEL = IER IER = LEVOLD IF (LEVEL.LT.0) LEVEL = 4 IF (LEVEL.GT.4) LEVEL = 4 GO TO 30 25 CONTINUE IF (LEVEL.LT.4) GO TO 30 C PRINT NON-DEFINED MESSAGE CALL UGETIO(1,NIN,IOUNIT) IF (IEQDF.EQ.1) WRITE(IOUNIT,50) IER,NAMEQ,IEQ,NAME IF (IEQDF.EQ.0) WRITE(IOUNIT,50) IER,NAME 30 IEQDF = 0 RETURN 35 FORMAT(19H *** TERMINAL ERROR,10X,7H(IER = ,I3, 1 20H) FROM IMSL ROUTINE ,A6,A1,A6) 40 FORMAT(36H *** WARNING WITH FIX ERROR (IER = ,I3, 1 20H) FROM IMSL ROUTINE ,A6,A1,A6) 45 FORMAT(18H *** WARNING ERROR,11X,7H(IER = ,I3, 1 20H) FROM IMSL ROUTINE ,A6,A1,A6) 50 FORMAT(20H *** UNDEFINED ERROR,9X,7H(IER = ,I5, 1 20H) FROM IMSL ROUTINE ,A6,A1,A6) C SAVE P FOR P = R CASE C P IS THE PAGE NAME C R IS THE ROUTINE NAME 55 IEQDF = 1 NAMEQ = NAME 65 RETURN END *DECK UGETIO C IMSL ROUTINE NAME - UGETIO C C----------------------------------------------------------------------- C C COMPUTER - CDC/SINGLE C C LATEST REVISION - JUNE 1, 1981 C C PURPOSE - TO RETRIEVE CURRENT VALUES AND TO SET NEW C VALUES FOR INPUT AND OUTPUT UNIT C IDENTIFIERS. C C USAGE - CALL UGETIO(IOPT,NIN,NOUT) C C ARGUMENTS IOPT - OPTION PARAMETER. (INPUT) C IF IOPT=1, THE CURRENT INPUT AND OUTPUT C UNIT IDENTIFIER VALUES ARE RETURNED IN NIN C AND NOUT, RESPECTIVELY. C IF IOPT=2, THE INTERNAL VALUE OF NIN IS C RESET FOR SUBSEQUENT USE. C IF IOPT=3, THE INTERNAL VALUE OF NOUT IS C RESET FOR SUBSEQUENT USE. C NIN - INPUT UNIT IDENTIFIER. C OUTPUT IF IOPT=1, INPUT IF IOPT=2. C NOUT - OUTPUT UNIT IDENTIFIER. C OUTPUT IF IOPT=1, INPUT IF IOPT=3. C C PRECISION/HARDWARE - SINGLE/ALL C C REQD. IMSL ROUTINES - NONE REQUIRED C C NOTATION - INFORMATION ON SPECIAL NOTATION AND C CONVENTIONS IS AVAILABLE IN THE MANUAL C INTRODUCTION OR THROUGH IMSL ROUTINE UHELP C C REMARKS EACH IMSL ROUTINE THAT PERFORMS INPUT AND/OR OUTPUT C OPERATIONS CALLS UGETIO TO OBTAIN THE CURRENT UNIT C IDENTIFIER VALUES. IF UGETIO IS CALLED WITH IOPT=2 OR C IOPT=3, NEW UNIT IDENTIFIER VALUES ARE ESTABLISHED. C SUBSEQUENT INPUT/OUTPUT IS PERFORMED ON THE NEW UNITS. C C COPYRIGHT - 1978 BY IMSL, INC. ALL RIGHTS RESERVED. C C WARRANTY - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN C APPLIED TO THIS CODE. NO OTHER WARRANTY, C EXPRESSED OR IMPLIED, IS APPLICABLE. C C----------------------------------------------------------------------- C SUBROUTINE UGETIO(IOPT,NIN,NOUT) C SPECIFICATIONS FOR ARGUMENTS INTEGER IOPT,NIN,NOUT C SPECIFICATIONS FOR LOCAL VARIABLES INTEGER NIND,NOUTD C FIRST EXECUTABLE STATEMENT IF (IOPT.EQ.3) GO TO 10 IF (IOPT.EQ.2) GO TO 5 IF (IOPT.NE.1) GO TO 9005 NIN = I1MACH(1) NOUT = I1MACH(2) GO TO 9005 5 NIND = NIN GO TO 9005 10 NOUTD = NOUT 9005 RETURN END *DECK MERFI C IMSL ROUTINE NAME - MERFI C C----------------------------------------------------------------------- C C COMPUTER - CDC/SINGLE C C LATEST REVISION - JANUARY 1, 1978 C C PURPOSE - INVERSE ERROR FUNCTION C C USAGE - CALL MERFI (P,Y,IER) C C ARGUMENTS P - INPUT VALUE IN THE EXCLUSIVE RANGE (-1.0,1.0) C Y - OUTPUT VALUE OF THE INVERSE ERROR FUNCTION C IER - ERROR PARAMETER (OUTPUT) C TERMINAL ERROR C IER = 129 INDICATES P LIES OUTSIDE THE LEGAL C RANGE. PLUS OR MINUS MACHINE INFINITY IS C GIVEN AS THE RESULT (SIGN IS THE SIGN OF C THE FUNCTION VALUE OF THE NEAREST LEGAL C ARGUMENT). C C PRECISION/HARDWARE - SINGLE/ALL C C REQD. IMSL ROUTINES - UERTST,UGETIO C C NOTATION - INFORMATION ON SPECIAL NOTATION AND C CONVENTIONS IS AVAILABLE IN THE MANUAL C INTRODUCTION OR THROUGH IMSL ROUTINE UHELP C C COPYRIGHT - 1978 BY IMSL, INC. ALL RIGHTS RESERVED. C C WARRANTY - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN C APPLIED TO THIS CODE. NO OTHER WARRANTY, C EXPRESSED OR IMPLIED, IS APPLICABLE. C C----------------------------------------------------------------------- C SUBROUTINE MERFI (P,Y,IER) C SPECIFICATIONS FOR ARGUMENTS REAL P,Y INTEGER IER C SPECIFICATIONS FOR LOCAL VARIABLES REAL A(65),H1,H2,H3,H4,X,SIGMA,Z,W,X3,X4,X5, * X6,B INTEGER N,IPP,L,LB2 DATA A(1),A(2),A(3),A(4),A(5),A(6),A(7),A(8),A(9), * A(10),A(11),A(12),A(13),A(14),A(15),A(16), * A(17),A(18),A(19),A(20),A(21),A(22),A(23) * /.99288537662,.12046751614, * .16078199342E-01,.26867044372E-02, * .49963473024E-03,.98898218599E-04, * .20391812764E-04,.4327271618E-05, * .938081413E-06,.206734721E-06, * .46159699E-07,.10416680E-07, * .2371501E-08,.543928E-09, * .125549E-09,.29138E-10, * .6795E-11,.1591E-11, * .374E-12,.88E-13, * .21E-13,.5E-14, * .1E-14/ DATA A(24),A(25),A(26),A(27),A(28),A(29),A(30), * A(31),A(32),A(33),A(34),A(35),A(36),A(37), * A(38),A(39),A(40) * /.91215880342E00,-.16266281868E-01, * .43355647295E-03,.21443857007E-03, * .2625751076E-05,-.302109105E-05, * -.12406062E-07,.62406609E-07, * -.540125E-09,-.142328E-08, * .34384E-10,.33585E-10, * -.1458E-11,-.81E-12, * .53E-13,.2E-13, * -.2E-14/ DATA A(41),A(42),A(43),A(44),A(45),A(46),A(47), * A(48),A(49),A(50),A(51),A(52),A(53),A(54), * A(55),A(56),A(57),A(58),A(59),A(60),A(61), * A(62),A(63),A(64),A(65) * /.95667970902,-.023107004309, * -.43742360975E-02,-.57650342265E-03, * -.10961022307E-04,.25108547025E-04, * .10562336068E-04,.275441233E-05, * .432484498E-06,-.20530337E-07, * -.43891537E-07,-.1768401E-07, * -.3991289E-08,-.186932E-09, * .272923E-09,.132817E-09, * .31834E-10,.1670E-11, * -.2036E-11,-.965E-12, * -.22E-12,-.1E-13, * .13E-13,.6E-14, * .1E-14/ DATA H1,H2,H3,H4/-1.5488130424, * 2.5654901231,-.55945763133, * 2.2879157163/ C FIRST EXECUTABLE STATEMENT X = P IER = 0 SIGMA = SIGN(1.,X) IF (.NOT.(X.GT.-1..AND.X.LT.1.)) GO TO 35 Z = ABS(X) IF(Z.GT..8) GO TO 20 W = Z*Z/.32-1. N = 22 IPP = 1 L = 1 5 LB2 = 1 X3 = 1. X4 = W X6 = A(IPP) 10 X6 = X6 + A(IPP+LB2) * X4 X5 = X4 * W * 2.-X3 X3 = X4 X4 = X5 LB2 = LB2 + 1 IF (LB2 .LE. N) GO TO 10 GO TO (15,30),L 15 Y = Z * X6 * SIGMA GO TO 9005 20 B = SQRT(-ALOG(1.-Z*Z)) IF (Z .GT. .9975) GO TO 25 W = H1*B+H2 IPP = 24 L = 2 N = 16 GO TO 5 25 W = H3 * B + H4 IPP = 41 N = 24 L = 2 GO TO 5 30 Y = B * X6 * SIGMA GO TO 9005 35 Y = SIGMA*R1MACH(2) IER = 129 9000 CONTINUE CALL UERTST(IER,'MERFI ') 9005 RETURN END *DECK MDNRIS C IMSL ROUTINE NAME - MDNRIS C C----------------------------------------------------------------------- C C COMPUTER - CDC/SINGLE C C LATEST REVISION - JUNE 1, 1981 C C PURPOSE - INVERSE STANDARD NORMAL (GAUSSIAN) C PROBABILITY DISTRIBUTION FUNCTION C C USAGE - CALL MDNRIS (P,Y,IER) C C ARGUMENTS P - INPUT VALUE IN THE EXCLUSIVE RANGE (0.0,1.0) C Y - OUTPUT VALUE OF THE INVERSE NORMAL (0,1) C PROBABILITY DISTRIBUTION FUNCTION C IER - ERROR PARAMETER (OUTPUT) C TERMINAL ERROR C IER = 129 INDICATES P LIES OUTSIDE THE LEGAL C RANGE. PLUS OR MINUS MACHINE INFINITY IS C GIVEN AS THE RESULT (SIGN IS THE SIGN OF C THE FUNCTION VALUE OF THE NEAREST LEGAL C ARGUMENT). C C PRECISION/HARDWARE - SINGLE/ALL C C REQD. IMSL ROUTINES - MERFI,UERTST,UGETIO C C NOTATION - INFORMATION ON SPECIAL NOTATION AND C CONVENTIONS IS AVAILABLE IN THE MANUAL C INTRODUCTION OR THROUGH IMSL ROUTINE UHELP C C COPYRIGHT - 1978 BY IMSL, INC. ALL RIGHTS RESERVED. C C WARRANTY - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN C APPLIED TO THIS CODE. NO OTHER WARRANTY, C EXPRESSED OR IMPLIED, IS APPLICABLE. C C----------------------------------------------------------------------- C SUBROUTINE MDNRIS (P,Y,IER) C SPECIFICATIONS FOR ARGUMENTS REAL P,Y INTEGER IER C SPECIFICATIONS FOR LOCAL VARIABLES REAL D(25),A,B,W,H3,H4,X3,X4,X5,X6 REAL SIGMA,SQRT2,X DATA SQRT2/1.4142135623731/ DATA D(1)/.95667970902049/ DATA D(2)/-.02310700430907/ DATA D(3)/-.00437423609751/ DATA D(4)/-.00057650342265/ DATA D(5)/-.00001096102231/ DATA D(6)/.00002510854703/ DATA D(7)/.00001056233607/ DATA D(8)/.00000275441233/ DATA D(9)/.00000043248450/ DATA D(10)/-.00000002053034/ DATA D(11)/-.00000004389154/ DATA D(12)/-.00000001768401/ DATA D(13)/-.00000000399129/ DATA D(14)/-.00000000018693/ DATA D(15)/.00000000027292/ DATA D(16)/.00000000013282/ DATA D(17)/.00000000003183/ DATA D(18)/.00000000000167/ DATA D(19)/-.00000000000204/ DATA D(20)/-.00000000000097/ DATA D(21)/-.00000000000022/ DATA D(22)/-.00000000000001/ DATA D(23)/.00000000000001/ DATA D(24)/.00000000000001/ DATA D(25)/.00000000000000/ DATA H3/-.55945763132983/ DATA H4/2.2879157162634/ C FIRST EXECUTABLE STATEMENT IER = 0 IF (P .GT. 0.0 .AND. P .LT. 1.0) GO TO 5 IER = 129 SIGMA = SIGN(1.0,P) Y = SIGMA * R1MACH(2) GO TO 9000 5 IF(P.LE.R1MACH(4)) GO TO 10 X = 1.0 -(P + P) CALL MERFI (X,Y,IER) Y = -SQRT2 * Y GO TO 9005 C P TOO SMALL, COMPUTE Y DIRECTLY 10 A = P+P B = SQRT(-ALOG(A+(A-A*A))) W = H3 * B + H4 X3 = 1.0 X4 = W X6 = D(1) DO 15 I=2,25 X6 = X6 + D(I) * X4 X5 = X4 * W * 2. - X3 X3 = X4 X4 = X5 15 CONTINUE Y = B*X6 Y = -SQRT2 * Y GO TO 9005 9000 CONTINUE CALL UERTST(IER,'MDNRIS') 9005 RETURN END *DECK MDCH C IMSL ROUTINE NAME - MDCH C C----------------------------------------------------------------------- C C COMPUTER - CDC/SINGLE C C LATEST REVISION - JUNE 1, 1980 C C PURPOSE - CHI-SQUARED PROBABILITY DISTRIBUTION FUNCTION C C USAGE - CALL MDCH (CS,DF,P,IER) C C ARGUMENTS CS - INPUT VALUE FOR WHICH THE PROBABILITY IS C COMPUTED. CS MUST BE GREATER THAN OR EQUAL C TO ZERO. C DF - INPUT NUMBER OF DEGREES OF FREEDOM OF THE C CHI-SQUARED DISTRIBUTION. DF MUST BE GREATER C THAN OR EQUAL TO .5 AND LESS THAN OR EQUAL C TO 200,000. C P - OUTPUT PROBABILITY THAT A RANDOM VARIABLE C WHICH FOLLOWS THE CHI-SQUARED DISTRIBUTION C WITH DF DEGREES OF FREEDOM IS LESS THAN OR C EQUAL TO CS. C IER - ERROR PARAMETER. (OUTPUT) C TERMINAL ERROR C IER = 129 INDICATES THAT CS OR DF WAS C SPECIFIED INCORRECTLY. C WARNING ERROR C IER = 34 INDICATES THAT THE NORMAL PDF C WOULD HAVE PRODUCED AN UNDERFLOW. C C PRECISION/HARDWARE - SINGLE/ALL C C REQD. IMSL ROUTINES - H32/MDNOR,MERRC=ERFC,MGAMAD=DGAMMA,UERTST, C UGETIO C - H36,H48,H60/MDNOR,MERRC=ERFC,MGAMA=GAMMA, C UERTST,UGETIO C C NOTATION - INFORMATION ON SPECIAL NOTATION AND C CONVENTIONS IS AVAILABLE IN THE MANUAL C INTRODUCTION OR THROUGH IMSL ROUTINE UHELP C C COPYRIGHT - 1978 BY IMSL, INC. ALL RIGHTS RESERVED. C C WARRANTY - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN C APPLIED TO THIS CODE. NO OTHER WARRANTY, C EXPRESSED OR IMPLIED, IS APPLICABLE. C C----------------------------------------------------------------------- C SUBROUTINE MDCH (CS,DF,P,IER) C SPECIFICATIONS FOR ARGUMENTS REAL CS,DF,P C SPECIFICATIONS FOR LOCAL VARIABLES REAL A,Z,DGAM,EPS,W,W1,B,Z1,ONE,HALF,X,C REAL GAMMA,THRD,PT2 DATA EPS/1.E-14/,HALF/5.E-1/ DATA THRTEN/13.E0/,ONE/1.E0/ DATA THRD/.33333333333333E0/ DATA PT2/.22222222222222E0/ FUNC(W,A,Z)=W* EXP(A*ALOG(Z)-Z) C FIRST EXECUTABLE STATEMENT C TEST FOR INVALID INPUT VALUES IF (DF .GE. .5 .AND. DF .LE. 2.E5 .AND. CS .GE. 0.0) GO TO 5 IER=129 P=-R1MACH(1) GO TO 9000 5 IER=0 C SET P=0. IF CS IS LESS THAN OR C EQUAL TO 10.**(-12) IF (CS .GT. 1.E-12) GO TO 15 10 P=0.0 GO TO 9005 15 IF(DF.LE.100.) GO TO 20 C USE NORMAL DISTRIBUTION APPROXIMATION C FOR LARGE DEGREES OF FREEDOM IF(CS.LT.2.0) GO TO 10 X=((CS/DF)**THRD-(ONE-PT2/DF))/SQRT(PT2/DF) IF (X .GT. 5.0) GO TO 50 IF (X .LT. -11.313708) GO TO 55 CALL MDNOR (X,P) GO TO 9005 C INITIALIZATION FOR CALCULATION USING C INCOMPLETE GAMMA FUNCTION 20 IF (CS .GT. 200.) GO TO 50 A=HALF*DF Z=HALF*CS C = A DGAM = GAMMA(C) W=AMAX1(HALF*A,THRTEN) IF (Z .GE. W) GO TO 35 IF (DF .GT. 25. .AND. CS .LT. 2.) GO TO 10 C CALCULATE USING EQUATION NO. 6.5.29 W=ONE/(DGAM*A) W1=W DO 25 I=1,50 B=I W1=W1*Z/(A+B) IF (W1 .LE. EPS*W) GO TO 30 W=W+W1 25 CONTINUE 30 P=FUNC(W,A,Z) GO TO 9005 C CALCULATE USING EQUATION NO. 6.5.32 35 Z1=ONE/Z B=A-ONE W1=B*Z1 W=ONE+W1 DO 40 I=2,50 B=B-ONE W1=W1*B*Z1 IF (W1 .LE. EPS*W) GO TO 45 W=W+W1 40 CONTINUE 45 W=Z1*FUNC(W,A,Z) P=ONE-W/DGAM GO TO 9005 50 P=1.0 GO TO 9005 C WARNING ERROR - UNDERFLOW WOULD HAVE C OCCURRED 55 P=0.0 IER=34 9000 CONTINUE CALL UERTST (IER,'MDCH ') 9005 RETURN END integer function i1mach(i) c i/o unit numbers. c c i1mach( 1) = the standard input unit. c c i1mach( 2) = the standard output unit. c c i1mach( 3) = the standard punch unit. c c i1mach( 4) = the standard error message unit. c c words. c c i1mach( 5) = the number of bits per word. c c i1mach( 6) = the number of characters per word. c c integers. c c assume integers are represented in the s-digit, base-a form c c sign ( x(s-1)*a**(s-1) + ... + x(1)*a + x(0) ) c c where 0 .le. x(i) .lt. a for i=0,...,s-1. c c i1mach( 7) = a, the base. c c i1mach( 8) = s, the number of base-a digits. c c i1mach( 9) = a**s - 1, the largest magnitude. c c floating-point numbers. c c assume floating-point numbers are represented in the t-digit, c base-b form c c sign (b**e)*( (x(1)/b) + ... + (x(t)/b**t) ) c c where 0 .lt. x(i) .lt. b for i=1,...,t, c 0 .lt. x(1), and emin .le. e .le. emax. c c i1mach(10) = b, the base. c c single-precision c c i1mach(11) = t, the number of base-b digits. c c i1mach(12) = emin, the smallest exponent e. c c i1mach(13) = emax, the largest exponent e. c c double-precision c c i1mach(14) = t, the number of base-b digits. c c i1mach(15) = emin, the smallest exponent e. c c i1mach(16) = emax, the largest exponent e. c c to alter this function for a particular environment, c the desired set of data statements should be activated by c removing the c from column 1. Also, the values of c i1mach(1) - i1mach(4) should be checked for consistency c with the local operating system. c integer imach(16),output c equivalence (imach(4),output) c c machine constants for the honeywell 6000 series. c c data imach( 1) / 5/ c data imach( 2) / 6/ c data imach( 3) / 43/ c data imach( 4) / 6/ c data imach( 5) / 36/ c data imach( 6) / 6/ c data imach( 7) / 2/ c data imach( 8) / 35/ c data imach( 9) /O377777777777/ c data imach(10) / 2/ c data imach(11) / 27/ c data imach(12) /-127/ c data imach(13) / 127/ c data imach(14) / 63/ c data imach(15) /-127/ c data imach(16) / 127/ c c machine constants for the ibm 360 and 370 series, c and the sel systems 85/86. c c data imach( 1) / 5/ c data imach( 2) / 6/ c data imach( 3) / 7/ c data imach( 4) / 6/ c data imach( 5) / 32/ c data imach( 6) / 4/ c data imach( 7) / 2/ c data imach( 8) / 31/ c data imach( 9) /Z7FFFFFFF/ c data imach(10) / 16/ c data imach(11) / 6/ c data imach(12) / -64/ c data imach(13) / 63/ c data imach(14) / 14/ c data imach(15) / -64/ c data imach(16) / 63/ c c machine constants for the cray 1. c c data imach( 1) /5linput / c data imach( 2) /6loutput/ c data imach( 3) /5lpunch / c data imach( 4) /6loutput/ c data imach( 5) / 64 / c data imach( 6) / 8 / c data imach( 7) / 2/ c data imach( 8) / 23/ c data imach( 9) /37777777B / c data imach(10) / 2 / c data imach(11) / 47 / c data imach(12) /-8192 / c data imach(13) / 8191 / c data imach(14) / 95 / c data imach(15) /-8144 / c data imach(16) / 8191 / c c machine constants for the cdc 6000 and 7000 series. c c for imach(2) on a cdc 7600 running ltss, c the initial value should be 3ltty. c c data imach( 1) / 5linput / c data imach( 2) / 6loutput/ c data imach( 3) / 5lpunch / c data imach( 4) / 6loutput / c data imach( 5) / 60/ c data imach( 6) / 10/ c data imach( 7) / 2/ c data imach( 8) / 48/ c data imach( 9) /00007777777777777777B/ c data imach(10) / 2/ c data imach(11) / 48/ c data imach(12) / -974/ c data imach(13) / 1070/ c data imach(14) / 96/ c data imach(15) / -926/ c data imach(16) / 1070/ c c machine constants for the pdp-10 (ka processor). c c data imach( 1) / 5/ c data imach( 2) / 6/ c data imach( 3) / 5/ c data imach( 4) / 6/ c data imach( 5) / 36/ c data imach( 6) / 5/ c data imach( 7) / 2/ c data imach( 8) / 35/ c data imach( 9) / ;377777777777/ c data imach(10) / 2/ c data imach(11) / 27/ c data imach(12) /-128/ c data imach(13) / 127/ c data imach(14) / 54/ c data imach(15) /-128/ c data imach(16)/ 127/ c c machine constants for the pdp-10 (ki processor). c c data imach( 1) / 5/ c data imach( 2) / 6/ c data imach( 3) / 5/ c data imach( 4) / 6/ c data imach( 5) / 36/ c data imach( 6) / 5/ c data imach( 7) / 2/ c data imach( 8) / 35/ c data imach( 9) / ;377777777777/ c data imach(10) / 2/ c data imach(11) / 27/ c data imach(12) /-128/ c data imach(13) / 127/ c data imach(14) / 62/ c data imach(15) /-128/ c data imach(16) / 127/ c c machine constants for the pdp-11. c c data imach( 1) / 5/ c data imach( 2) / 6/ c data imach( 3) / 5/ c data imach( 4) / 6/ c data imach( 5) / 32/ c data imach( 6) / 4/ c data imach( 7) / 2/ c data imach( 8) / 31/ c data imach( 9) / 2147483647/ c data imach(10) / 2/ c data imach(11) / 24/ c data imach(12) / -127/ c data imach(13) / 127/ c data imach(14) / 56/ c data imach(15) / -127/ c data imach(16) / 127/ c c machine constants for the univac 1100 series. c c data imach( 1) / 5/ c data imach( 2) / 6/ c data imach( 3) / 7/ c data imach( 4) / 6/ c data imach( 5) / 36/ c data imach( 6) / 6/ c data imach( 7) / 2/ c data imach( 8) / 35/ c data imach( 9) /O377777777777/ c data imach(10) / 2/ c data imach(11) / 27/ c data imach(12) /-128/ c data imach(13) / 127/ c data imach(14) / 61/ c data imach(15) /-1024/ c data imach(16) / 1023/ c c c machine constants for IBM PC Professional Fortran Compiler c data imach( 1) / 5/ data imach( 2) / 6/ data imach( 3) / 7/ data imach( 4) / 0/ data imach( 5) / 32/ data imach( 6) / 4/ data imach( 7) / 2/ data imach( 8) / 31/ data imach( 9) /2147483647/ data imach(10) / 2/ data imach(11) / 23/ data imach(12) /-125/ data imach(13) / 127/ data imach(14) / 52/ data imach(15) /-1021/ data imach(16) /1023/ c if (i .lt. 1 .or. i .gt. 16) go to 10 c i1mach=imach(i) return c 10 write(output,9000) 9000 format('1error 1 in i1mach - i out of bounds') c c call fdump c stop 'I1MACH bombed' c end c real function r1mach(i) c c single-precision machine constants c c r1mach( 1) = b**(emin-1), the smallest positive magnitude. c c r1mach( 2) = b**emax*(1 - b**(-t)), the largest magnitude. c c r1mach( 3) = b**(-t), the smallest relative spacing. c c r1mach( 4) = b**(1-t), the largest relative spacing. c c r1mach( 5) = log10(b) c c to alter this function for a particular environment, c the desired set of data statements should be activated by c removing the c from column 1. c c where possible, octal or hexadecimal constants have been used c to specify the constants exactly. c real rmach(5) integer output c c machine constants for the honeywell 6000 series. c c data rmach( 1) / O402400000000/ c data rmach( 2) / O376777777777/ c data rmach( 3) / O714400000000/ c data rmach( 4) / O716400000000/ c data rmach( 5) / O776464202324/ c c machine constants for the IBM 360 and 370 series, c and the sel systems 85/86. c c data rmach( 1) /Z00100000/ c data rmach( 2) /Z7FFFFFFF/ c data rmach( 3) /Z3B100000/ c data rmach( 4) /Z3C100000/ c data rmach( 5) /Z41134413/ c c machine constants for the Cray 1. c c data rmach( 1) / 200004000000000000000B / c data rmach( 2) / 577777777777777777777B / c data rmach( 3) / 377224000000000000000B / c data rmach( 4) / 377234000000000000000B / c data rmach( 5) / 0.3010299956 63981195E0 / c rmach(3) = 2**(-47) = 7.10543E-15 c rmach(4) = 2**(-46) = 1.42109E-14 c c machine constants for the CDC 6000 and 7000 series. c c data rmach( 1) / 00014000000000000000B/ c data rmach( 2) / 37767777777777777777B/ c data rmach( 3) / 16404000000000000000B/ c data rmach( 4) / 16414000000000000000B/ c data rmach( 5) / 17164642023241175720B/ c c machine constants for the PDP-10 (KA or KI processor). c c data rmach( 1) / ;000400000000/ c data rmach( 2) / ;377777777777/ c data rmach( 3) / ;146400000000/ c data rmach( 4) / ;147400000000/ c data rmach( 5) / ;177464202324/ c c machine constants for the PDP-11. c c data rmach( 1) / 0.29387358771E-38/ c data rmach( 2) / 0.17014117331E+39/ c data rmach( 3) / 0.59604644775E-07/ c data rmach( 4) / 0.11920928955E-06/ c data rmach( 5) / 0.30102999566 / c c machine constants for the UNIVAC 1100 series. c c data rmach( 1) / O000400000000/ c data rmach( 2) / O377777777777/ c data rmach( 3) / O146400000000/ c data rmach( 4) / O147400000000/ c data rmach( 5) / O177464202324/ c c machine constants for the IBM PC Professional Fortran Compiler c data rmach( 1) / 1.1755E-38/ data rmach( 2) / 3.4028E+38/ data rmach( 3) / 1.1921E-07/ data rmach( 4) / 2.3842E-07/ data rmach( 5) / 3.0103E-01/ c if (i.lt.1 .or. i.gt.5) then output=i1mach(4) write(output,*) ' R1MACH - I out of bounds' stop 'R1MACH bombed' endif c r1mach = rmach(i) return c end double precision function d1mach(i) c c double-precision machine constants c c d1mach( 1) = b**(emin-1), the smallest positive magnitude. c c d1mach( 2) = b**emax*(1 - b**(-t)), the largest magnitude. c c d1mach( 3) = b**(-t), the smallest relative spacing. c c d1mach( 4) = b**(1-t), the largest relative spacing. c c d1mach( 5) = log10(b) c c to alter this function for a particular environment, c the desired set of data statements should be activated by c removing the c from column 1. c c where possible, octal or hexadecimal constants have been used c to specify the constants exactly which has in some cases c required the use of equivalent integer arrays. c integer small(2) integer large(2) integer right(2) integer diver(2) integer log10(2) c double precision dmach(5) integer output c equivalence (dmach(1),small(1)) equivalence (dmach(2),large(1)) equivalence (dmach(3),right(1)) equivalence (dmach(4),diver(1)) equivalence (dmach(5),log10(1)) c c machine constants for the honeywell 6000 series. c c data small(1),small(2) / O402400000000,O000000000000/ c data large(1),large(2) / O376777777777,O777777777777/ c data right(1),right(2) / O604400000000,O000000000000/ c data diver(1),diver(2) / O606400000000,O000000000000/ c data log10(1),log10(2) / O776464202324,O117571775714/ c c machine constants for the IBM 360 and 370 series, c and the sel systems 85/86. c c data small(1),small(2) /Z00100000,Z00000000/ c data large(1),large(2) /Z7FFFFFFF,ZFFFFFFFF/ c data right(1),right(2) /Z3B100000,Z00000000/ c data diver(1),diver(2) /Z3C100000,Z00000000/ c data log10(1),log10(2) /Z41134413,Z509F79FF/ c c machine constants for the Cray 1. c c data small(1),small(2) / 200604000000000000000B,20000000000000 c 10000000B/ c data large(1),large(2) / 577777777777777777777B,57717777777777 c 17777777B/ c data right(1),right(2) / 376424000000000000000B,37562000000000 c 10000000B/ c data diver(1),diver(2) / 376234000000000000000B,37563000000000 c 10000000B/ c data dmach( 5) / 0.3010299956 6398119521 3738894724D+00/ c note that dmach(3)=0.252435D-28 c and dmach(4)=0.5048709793D-28 c c machine constants for the CDC 6000 and 7000 series. c c data small(1),small(2) / 00014000000000000000B, c 1 00010000000000000000B/ c data large(1),large(2) / 37767777777777777777B, c 1 37167777777777777777B/ c data right(1),right(2) / 16404000000000000000B, c 1 15000000000000000000B/ c data diver(1),diver(2) / 16414000000000000000B, c 1 15010000000000000000B/ c data log10(1),log10(2) / 17164642023241175720B, c 1 16367571421742254654B/ c c machine constants for the PDP-10 (KA processor). c c data small(1),small(2) / ;000400000000,;000000000000/ c data large(1),large(2) / ;377777777777,;344777777777/ c data right(1),right(2) / ;146400000000,;000000000000/ c data diver(1),diver(2) / ;147400000000,;000000000000/ c data log10(1),log10(2) / ;177464202324,;144117571776/ c c machine constants for the PDP-10 (KA processor). c c data small(1),small(2) / ;000400000000,;000000000000/ c data large(1),large(2) / ;377777777777,;344777777777/ c data right(1),right(2) / ;103400000000,;000000000000/ c data diver(1),diver(2) / ;104400000000,;000000000000/ c data log10(1),log10(2) / ;177464202324,;476747767461/ c c machine constants for the PDP-11. c c data dmach( 1) / 0.2938735877055719 D-38/ c data dmach( 2) / 0.170141183460469229D+39/ c data dmach( 3) / 0.138777878078144568D-16/ c data dmach( 4) / 0.277555756156289135D-16/ c data dmach( 5) / 0.301029995663981195D+00/ c c machine constants for the UNIVAC 1100 series. c c data small(1),small(2) / O000400000000,O000000000000/ c data large(1),large(2) / O377777777777,O777777777777/ c data right(1),right(2) / O170540000000,O000000000000/ c data diver(1),diver(2) / O170640000000,O000000000000/ c data log10(1),log10(2) / O177746420232,O411757177572/ c c machine constants for the IBM PC Professional Fortran Compiler c data dmach( 1) / 2.23D-308/ data dmach( 2) / 1.79D+308/ data dmach( 3) / 2.2204D-52/ data dmach( 4) / 4.4409D-16/ data dmach( 5) / 3.0103D-01/ c if (i.lt.1 .or. i.gt.5) then output=i1mach(4) write(output,*) ' D1MACH - I out of bounds' stop 'D1MACH bombed' endif c d1mach = dmach(i) return c end