PROGRAM TSCALE C C ROBERT H. BLESSING C MEDICAL FOUNDATION OF BUFFALO C 73 HIGH STREET C BUFFALO, NEW YORK 14203, USA C TELEPHONE: (716) 856-9600 C ELECTRONIC MAIL: Blessing@MFB.Buffalo.Edu C C JULY 1993 C PARAMETER (MSTD=20,MCUT=10,MDAT=500,MOMIT=20,MREJ=50) CHARACTER TITLE*80,ATIME*8,ADATE*9,FILE1*40,FILE2*40 COMMON /BLOCK1/ TITLE,ATIME,ADATE COMMON /BLOCK2/ IO1,IO2,ILP,ITABL,IPLOT,IQSUM,IGRAPH, & IAPPLY,IANISO,IYDIFF,ISDIFF,A(6),GINV(3,3) COMMON /BLOCK3/ NSTD,IHKL(MSTD,3),NCUT(MSTD),XCUT(MSTD,MCUT),XZERO COMMON /BLOCK4/ NIOMIT(MSTD),IOMIT(MSTD,MOMIT,2), & NXOMIT(MSTD),XOMIT(MSTD,MOMIT,2) COMMON /BLOCK5/ NIREJ,IREJ(MREJ,2), & NXREJ,XREJ(MREJ,2) COMMON /BLOCK6/ NDAT(MSTD,MCUT),XX(MSTD,MCUT,MDAT), & YY(MSTD,MCUT,MDAT),SS(MSTD,MCUT,MDAT) COMMON /BLOCK7/ PAR(MSTD,MCUT,4),COV(MSTD,MCUT,4,4), & XLIM(MSTD,MCUT,2),YZERO(MSTD),P CALL INPUT CALL CULL CALL ANALYZE CALL APPLY STOP 'PROGRAM TSCALE FINIS!' END C----------------------------------------------------------------------- SUBROUTINE ANALYZE PARAMETER (MSTD=20,MCUT=10,MDAT=500,MOMIT=20,MREJ=50) PARAMETER (NX=101,NY=51) CHARACTER TITLE*80,ATIME*8,ADATE*9 COMMON /BLOCK1/ TITLE,ATIME,ADATE COMMON /BLOCK2/ IO1,IO2,ILP,ITABL,IPLOT,IQSUM,IGRAPH, & IAPPLY,IANISO,IYDIFF,ISDIFF,A(6),GINV(3,3) COMMON /BLOCK3/ NSTD,IHKL(MSTD,3),NCUT(MSTD),XCUT(MSTD,MCUT),XZERO COMMON /BLOCK6/ NDAT(MSTD,MCUT),XX(MSTD,MCUT,MDAT), & YY(MSTD,MCUT,MDAT),SS(MSTD,MCUT,MDAT) COMMON /BLOCK7/ PAR(MSTD,MCUT,4),COV(MSTD,MCUT,4,4), & XLIM(MSTD,MCUT,2),YZERO(MSTD),P DIMENSION X(MDAT),Y(MDAT),S(MDAT),C(4),CC(4,4) DIMENSION SZERO(MSTD),PP(MSTD,MCUT) C C ANALYZE STANDARD REFLECTION INTENSITIES. C DO 10 I=1,NSTD DO 10 J=1,NCUT(I)+1 N=NDAT(I,J) IF (N.EQ.0) GO TO 10 DO 50 K=1,4 PAR(I,J,K)=0 DO 50 L=1,4 COV(I,J,K,L)=0 50 CONTINUE XMIN=+1E6 XMAX=-1E6 YMIN=+1E6 YMAX=-1E6 SUMN=0 SUMW=0 SUMWY=0 SUMWY2=0 C C PICK OUT DATA AND EVALUATE STATISTICS FOR JTH SEGMENT OF ITH C STANDARD. C DO 11 K=1,N X(K)=XX(I,J,K) Y(K)=YY(I,J,K) S(K)=SS(I,J,K) XI=X(K) YI=Y(K) SIGYI=S(K) WI=1/SIGYI**2 SUMN=SUMN+1 SUMW=SUMW+WI SUMWY=SUMWY+WI*YI SUMWY2=SUMWY2+WI*YI**2 XMIN=MIN(XMIN,XI) XMAX=MAX(XMAX,XI) 11 CONTINUE XLIM(I,J,1)=XMIN XLIM(I,J,2)=XMAX YMEAN=SUMWY/SUMW RMSD=SQRT(ABS(SUMWY2/SUMW-YMEAN**2)) ESD=SQRT(1/SUMW) IF (SUMN.LT.2) GO TO 12 RMSD=RMSD*SQRT(SUMN/(SUMN-1)) ESD=ESD*SQRT(SUMN) 12 IF (ITABL.EQ.0) GO TO 13 WRITE (ILP,102) ATIME,ADATE,TITLE WRITE (ILP,101) (IHKL(I,K),K=1,3),YMEAN,ESD,RMSD, & (X(K),Y(K),K=1,N) 13 CONTINUE C C FIT A CURVE TO THE DATA (X,Y) FOR THE SEGMENT. C IF (N.GT.2) GO TO 14 C C FOR SEGMENTS WITH ONLY ONE OR TWO DATA POINTS, ASSUME Y(X) = C CONSTANT. C C(1)=YMEAN CC(1,1)=MAX(RMSD,ESD) DO 130 K=2,4 C(K)=0 DO 130 L=2,4 CC(K,L)=0 130 CONTINUE GO TO 140 14 CONTINUE C C FOR SEGMENTS WITH THREE OR MORE DATA POINTS, FIT A POWER C SERIES POLYNOMIAL, Y(X) = A0 + A1*X + A2*X**2 + A3*X**3, TO C THE SEGMENT. C CALL REGRESS (N,X,Y,S,C,CC) 140 DO 141 K=1,4 PAR(I,J,K)=C(K) DO 141 L=1,4 COV(I,J,K,L)=CC(K,L) 141 CONTINUE C C CALCULATE PROPORTIONALITY FACTOR, P, FOR INSTRUMENTAL C INSTABILITIES FOR THE SEGMENT. C CALL PCALC (N,X,Y,S,C,PSQ) PP(I,J)=PSQ C C PLOT DATA AND FITTED CURVE FOR THE SEGMENTS WITH TWO OR MORE C DATA POINTS. C IF (N.LT.2) GO TO 10 IF (IPLOT.EQ.0) GO TO 17 DO 15 K=1,N YI=Y(K) SIGYI=S(K) YMIN=MIN(YMIN,YI-2*SIGYI) YMAX=MAX(YMAX,YI+2*SIGYI) 15 CONTINUE YMID=(YMIN+YMAX)/2 DELY=(NY-1)*ESD Y1=YMID-DELY/2 Y2=YMID+DELY/2 YMIN=MIN(YMIN,Y1) YMAX=MAX(YMAX,Y2) DELX=(XMAX-XMIN)/(NX-1) DELY=(YMAX-YMIN)/(NY-1) WRITE (ILP,110) ATIME,ADATE,TITLE CALL PLOT1 (N,X,Y,XMIN,XMAX,YMIN,YMAX,C) WRITE (ILP,100) (IHKL(I,K),K=1,3),(PAR(I,J,K),K=1,4), & (SQRT(COV(I,J,K,K)),K=1,4),SQRT(MAX(PP(I,J),0.0)),DELX,DELY,ESD 17 CONTINUE C C EVALUATE AND PLOT Q-SUMS FOR THE SEGMENT. C IF (IQSUM.EQ.0) GO TO 10 QSUM=0 YMIN=+1E6 YMAX=-1E6 DO 16 K=1,N QSUM=QSUM+(Y(K)-YMEAN) Y(K)=QSUM YMIN=MIN(YMIN,QSUM) YMAX=MAX(YMAX,QSUM) 16 CONTINUE YMID=(YMIN+YMAX)/2 DELY=(NY-1)*ESD Y1=YMID-DELY/2 Y2=YMID+DELY/2 YMIN=MIN(YMIN,Y1) YMAX=MAX(YMAX,Y2) DELY=(YMAX-YMIN)/(NY-1) WRITE (ILP,110) ATIME,ADATE,TITLE CALL PLOT2 (N,X,Y,XMIN,XMAX,YMIN,YMAX) WRITE (ILP,103) (IHKL(I,K),K=1,3),YMEAN,ESD,DELX,DELY 10 CONTINUE C C CALCULATE Y0 = Y(X0) AND SIGMA(Y0) FOR EACH STANDARD. C WRITE (ILP,102) ATIME,ADATE,TITLE WRITE (ILP,104) XZERO DO 30 I=1,NSTD IF (NCUT(I).EQ.0.AND.NDAT(I,1).EQ.0) GO TO 30 C C FIND THE SEGMENT FOR THE ITH STANDARD THAT INCLUDES X = X0 IN C ITS DOMAIN. C N=0 DO 31 J=1,NCUT(I)+1 IF (NDAT(I,J).EQ.0) GO TO 31 N=N+1 X1=XLIM(I,J,1) X2=XLIM(I,J,2) IF (X1.LE.XZERO.AND.XZERO.LE.X2) GO TO 32 31 CONTINUE IF (N.EQ.0) GO TO 30 C C IF NO SEGMENT OF THE ITH STANDARD INCLUDES X0, THEN RESORT TO C THE SEGMENT NEAREST X0, AND EXTRAPOLATE TO X0. C DMIN=1E6 DO 33 J=1,NCUT(I)+1 X1=XLIM(I,J,1) X2=XLIM(I,J,2) D1=ABS(X1-XZERO) D2=ABS(X2-XZERO) D=MIN(D1,D2) IF (D.GE.DMIN) GO TO 33 DMIN=D JMIN=J 33 CONTINUE J=JMIN 32 CONTINUE C C CALCULATE Y0 AND SIGMA(Y0). C YZERO(I)=0 DO 36 K=1,4 YZERO(I)=YZERO(I)+PAR(I,J,K)*XZERO**(K-1) 36 CONTINUE SZERO(I)=0 DO 37 K=1,4 DO 37 L=1,4 SZERO(I)=SZERO(I)+XZERO**(K+L-2)*COV(I,J,K,L) 37 CONTINUE SZERO(I)=SQRT(ABS(SZERO(I))) WRITE (ILP,105) I,(IHKL(I,K),K=1,3),YZERO(I),SZERO(I) 30 CONTINUE C C CONVERT TO NORMALIZED INVERSE SCALING FACTORS Y/Y0 FOR EACH C SEGMENT. C WRITE (ILP,107) DO 20 I=1,NSTD DO 20 J=1,NCUT(I)+1 IF (NDAT(I,J).EQ.0) GO TO 20 DO 21 K=1,4 PAR(I,J,K)=PAR(I,J,K)/YZERO(I) DO 21 L=1,K T=(COV(I,J,K,L)+PAR(I,J,K)*PAR(I,J,L)*SZERO(I)**2)/YZERO(I)**2 COV(I,J,K,L)=T COV(I,J,L,K)=T 21 CONTINUE DO 22 K=1,NDAT(I,J) YI=YY(I,J,K)/YZERO(I) SI=YI**2*((SS(I,J,K)/YY(I,J,K))**2+(SZERO(I)/YZERO(I))**2) SI=SQRT(SI) YY(I,J,K)=YI SS(I,J,K)=SI 22 CONTINUE WRITE (ILP,108) I,J,(IHKL(I,K),K=1,3),NDAT(I,J),(XLIM(I,J,K),K=1, & 2),(PAR(I,J,K),K=1,4),(SQRT(COV(I,J,K,K)),K=1,4) 20 CONTINUE C C EVALUATE N-WEIGHTED, SQUARED-INTENSITY-WEIGHTED ROOT-MEAN- C SQUARE P-VALUE C WRITE (ILP,102) ATIME,ADATE,TITLE WRITE (ILP,109) PSQMIN=+1 PSQMAX=-1 SUMW=0 SUMWPSQ=0 DO 25 I=1,NSTD DO 25 J=1,NCUT(I)+1 N=NDAT(I,J) IF (N.EQ.0) GO TO 25 Y0=YZERO(I) PSQ=PP(I,J) WRITE (ILP,120) I,J,(IHKL(I,K),K=1,3),N,Y0,PSQ,SQRT(MAX(PSQ,0.0)) PSQMIN=MIN(PSQMIN,PSQ) PSQMAX=MAX(PSQMAX,PSQ) W=N*Y0**2 SUMW=SUMW+W SUMWPSQ=SUMWPSQ+W*PSQ 25 CONTINUE PSQ=SUMWPSQ/SUMW P=SQRT(MAX(PSQ,0.0)) WRITE (ILP,106) P,PSQ,PSQMIN,PSQMAX C C PLOT NORMALIZED DATA POINTS AND INVERSE SCALING FUNCTIONS FOR C EACH STANDARD AND COMPOSITE INVERSE SCALING FUNCTION AND ITS C FRACTIONAL ERROR. C IF (IGRAPH.NE.0) CALL GRAPH RETURN 110 FORMAT('1'/'0',10X,'PROGRAM TSCALE. ',A,1X,A,'. ',A) 100 FORMAT('0',10X,'REFLECTION INTENSITIES. Y = A0 + A1*X + A2*X**2', &' + A3*X**3, WHERE Y = I/LP AND X = T(HR).'/ &'0',10X,' H K K A0 A1 A2 A3 SIG(A0 &) SIG(A1) SIG(A2) SIG(A3) P'/ &' ',10X,3I3,4E10.3,4E9.2,2x,E10.3/ &'0',10X,'DELTA(X) = ',F10.3,' HR AND DELTA(Y) = ',E10.3,' PER GRID & POINT. SIGMA(Y) = ',E10.3) 101 FORMAT('0','STANDARD REFERENCE REFLECTION INTENSITY DATA Y = ', & 'I/LP AS A FUNCTION OF X = T(HR):'/ & '0','YMEAN = SUM(W*Y)/SUM(W)'/ & ' ','ESD = SQRT(N/SUM(W))'/ & ' ','RMSD = SQRT((N/(N-1))*(SUM(W*Y**2)/SUM(W) - ', & 'YMEAN**2))'/ & '0',' H K L YMEAN ESD RMSD'/ & ' ',3I3,3F10.2/'0',5(9X,'X',9X,'Y',5X)/ & ('0',5(F10.3,F10.2,5X))) 102 FORMAT('1'/'0PROGRAM TSCALE. ',A,1X,A,'. ',A) 103 FORMAT('0',10X,'Q-SUMS. Q(X(J)) = SUM(I=1,J)(Y(X(I)) - YMEAN) ', & 'AS A FUNCTION OF X = T(HR).'/ & '0',10X,' H K L YMEAN ESD'/ & ' ',10X,3I3,2F10.2/ & '0',10X,'DELTA(X) = ',F10.3,' HR AND DELTA(Q) = ',E10.3, & ' PER GRID POINT.') 109 FORMAT('0','PROPORTIONALITY FACTORS FOR INSTRUMENTAL INSTABILITIES & FOR EACH SEGMENT:'/ & '0','PSQ = (SUM(DELTA(Y)**2) - SUM(SIGMA(Y)**2))/SUM(Y**2)' &/ & '0','PSQ = (SUM((YI - Y(XI))**2) - SUM(SIGMA(YI)**2)/SUM(YI &**2))'/ & '0',' I J H K L N Y0 PSQ', &' P = SQRT(MAX(PSQ,0.0))'/) 120 FORMAT(' ',2I3,2X,3I3,2X,I5,F10.2,E12.3,F10.5) 106 FORMAT('0','N-WEIGHTED, SQUARED-INTENSITY-WEIGHTED ROOT-MEAN-SQUAR &E P-VALUE:'/ & '0','P = ',F8.5/ & '0','PSQ = ',E10.3/ & '0','PSQMIN = ',E10.3/ & ' ','PSQMAX = ',E10.3) 104 FORMAT('0','REFERENCE INTENSITY VALUES Y=I/LP CALCULATED AT THE ', & 'REFERENCE TIME OF THE EXPERIMENT X0 = ',F10.3,' HR'/ & '0',' I H K L Y0=Y(X0) SIGMA(Y0)'/) 105 FORMAT(' ',I3,5X,3I3,2F10.2) 107 FORMAT('0','NORMALIZED INVERSE SCALING FACTORS. Y/Y0 = A0 + A1*', &'X + A2*X**2 + A3*X**3, WHERE X = T(HR).'/ &'0',' I J H K L N XMIN XMAX A0 ', &' A1 A2 A3 SIG(A0) SIG(A1) SIG(A2) SIG(A3)'/) 108 FORMAT(' ',2I3,2X,3I3,2X,I6,2F8.2,4E11.3,4E9.2) END C----------------------------------------------------------------------- SUBROUTINE APPLY C C SCALES REFLECTION DATA ACCORDING TO THE SCALING FUNCTION C FROM SUBROUTINE ANALYZE. C PARAMETER (MSTD=20,MCUT=10,MDAT=500,MOMIT=20,MREJ=50) CHARACTER TITLE*80,ATIME*8,ADATE*9 COMMON /BLOCK1/ TITLE,ATIME,ADATE COMMON /BLOCK2/ IO1,IO2,ILP,ITABL,IPLOT,IQSUM,IGRAPH, & IAPPLY,IANISO,IYDIFF,ISDIFF,A(6),GINV(3,3) COMMON /BLOCK3/ NSTD,IHKL(MSTD,3),NCUT(MSTD),XCUT(MSTD,MCUT),XZERO COMMON /BLOCK5/ NIREJ,IREJ(MREJ,2), & NXREJ,XREJ(MREJ,2) COMMON /BLOCK6/ NDAT(MSTD,MCUT),XX(MSTD,MCUT,MDAT), & YY(MSTD,MCUT,MDAT),SS(MSTD,MCUT,MDAT) COMMON /BLOCK7/ PAR(MSTD,MCUT,4),COV(MSTD,MCUT,4,4), & XLIM(MSTD,MCUT,2),YZERO(MSTD),P DIMENSION U(3),V(3) IF (IAPPLY.EQ.0) RETURN C C TEST WHETHER SCALING IS TO BE WEIGHTED ANISOTROPICALLY OR BY C INTENSITY OR SIN(THETA)/LAMBDA DIFFERENCES. C IF (IANISO.EQ.0.AND.IYDIFF.EQ.0.AND.ISDIFF.EQ.0) GO TO 40 C C PRINT MESSAGES. C WRITE (ILP,300) ATIME,ADATE,TITLE IF (IANISO.NE.0.OR.ISDIFF.NE.0) WRITE (ILP,400) (A(I),I=1,6), & ((GINV(I,J),J=1,3),I=1,3) IF (IANISO.NE.0) C1=IANISO IF (IYDIFF.NE.0) C2=IYDIFF IF (ISDIFF.NE.0) C3=ISDIFF IF (C1.LT.0) C1=-1/C1 IF (C2.LT.0) C2=-1/C2 IF (C3.LT.0) C3=-1/C3 IF (IANISO.NE.0) WRITE (ILP,401) C1 IF (IYDIFF.NE.0) WRITE (ILP,402) C2 IF (ISDIFF.NE.0) WRITE (ILP,403) C3 40 CONTINUE C C INITIALIZE SOME VARIABLES. C WRITE (ILP,309) ATIME,ADATE,TITLE NSUM=0 SUMF=0 SUME=0 FMIN=+1E6 EMIN=+1E6 FMAX=-1E6 EMAX=-1E6 NREAD=0 NWRIT=0 NPRINT=0 C C LOOP THROUGH REFLECTION DATA FILE. C REWIND IO1 IEND=0 1 CALL READ1 (IO1,IEND,II,IH,IK,IL,A1,A2,A3,A4,Y,SIGY,XTIME) IF (IEND.NE.0) GO TO 9 NREAD=NREAD+1 C C X**N IS UNDEFINED FOR X = N = 0. C IF (XTIME.EQ.0) XTIME=0.0001 X=XTIME C C IS THIS REFLECTION TO BE REJECTED FROM THE DATA SET? C IF (NIREJ.EQ.0) GO TO 25 DO 20 I=1,NIREJ I1=IREJ(I,1) I2=IREJ(I,2) IF (I1.LE.ABS(II).AND.ABS(II).LE.I2) GO TO 1 20 CONTINUE 25 IF (NXREJ.EQ.0) GO TO 29 DO 21 I=1,NXREJ X1=XREJ(I,1) X2=XREJ(I,2) IF (X1.LT.X.AND.X.LT.X2) GO TO 1 21 CONTINUE 29 CONTINUE C C FIND THE INVERSE SCALING FACTOR APPLICABLE AT XTIME = X BY C AVERAGING OVER THE POLYNOMIAL CURVE SEGMENTS WHOSE DOMAIN C INCLUDES X = XTIME. C SUMN=0 SUMW=0 SUMWF=0 SUMWF2=0 SUMW2V=0 C C Q IS AN EXPANSION FACTOR USED TO DEAL WITH THE (RARE) CASE OF C A MEASUREMENT TIME NOT IN THE DOMAIN OF ANY SEGMENT. C EXTRAPOLATE NO MORE THAN HALF THE AVERAGE TIME BETWEEN C MEASUREMENTS OF STANDARDS. C Q=0 10 IF (Q.GT.0.5) GO TO 1 DO 11 I=1,NSTD DO 11 J=1,NCUT(I)+1 IF (NDAT(I,J).EQ.0) GO TO 11 N=NDAT(I,J) X1=XLIM(I,J,1) X2=XLIM(I,J,2) DX=Q*(X2-X1)/N X1=X1-DX X2=X2+DX IF (X.LT.X1) GO TO 11 IF (X.GT.X2) GO TO 11 F=0 VAR=0 DO 12 K=1,4 F=F+PAR(I,J,K)*X**(K-1) DO 12 L=1,4 VAR=VAR+X**(K+L-2)*COV(I,J,K,L) 12 CONTINUE C C VARIANCE CAN BE ZERO OR, DUE TO NUMERICAL IMPRECISION IN C COV(L,M), HAVE A SMALL NEGATIVE VALUE WHEN X IS NEAR X = XZERO C VAR=MAX(VAR,1E-9) C C RELATIVE WEIGHT FOR AVERAGING W = W0*W1*W2*W3. C C W0 FROM VARIANCE DUE TO POLYNOMIAL COEFFICIENTS. C C W1 FROM ANGLE BETWEEN RECIPROCAL LATTICE VECTORS FOR GIVEN C REFLECTION AND REFERENCE REFLECTION. C C W2 FROM DIFFERENCE BETWEEN INTENSITIES OF GIVEN REFLECTION C AND REFERENCE REFLECTION. C C W3 FROM DIFFERENCE BETWEEN SIN(THETA)/LAMBDA VALUES OF GIVEN C REFLECTION AND REFERENCE REFLECTION. C W0=1/VAR W1=1 W2=1 W3=1 IF (IANISO.EQ.0.AND.ISDIFF.EQ.0) GO TO 50 U(1)=IHKL(I,1) U(2)=IHKL(I,2) U(3)=IHKL(I,3) V(1)=IH V(2)=IK V(3)=IL 50 CONTINUE C C FUNCTION VTMV CALCULATES THE BILINEAR FORM, C SCALAR = (VECTOR TRANSPOSE)*(MATRIX)*(VECTOR). C IF (IANISO.EQ.0) GO TO 51 COSPH=VTMV(U,GINV,V)/SQRT(VTMV(U,GINV,U)*VTMV(V,GINV,V)) W1=1/(1+C1**2*(1-COSPH**2)) 51 IF (IYDIFF.EQ.0) GO TO 52 DELTA=Y-YZERO(I)*F W2=1/(1+(C2*DELTA)**2) 52 IF (ISDIFF.EQ.0) GO TO 53 DELTA=SQRT(VTMV(U,GINV,U))-SQRT(VTMV(V,GINV,V)) W3=1/(1+(C3*DELTA)**2) 53 W=W0*W1*W2*W3 SUMN=SUMN+1 SUMW=SUMW+W SUMWF=SUMWF+W*F SUMWF2=SUMWF2+W*F**2 SUMW2V=SUMW2V+W**2*VAR 11 CONTINUE Q=Q+0.5 IF (SUMN.EQ.0) GO TO 10 C C CALCULATE MEAN INVERSE SCALING FACTOR AND ITS STANDARD C DEVIATION. C F=SUMWF/SUMW WMSD=SUMWF2/SUMW-F**2 IF (SUMN.GT.1) WMSD=WMSD/(SUMN-1) VAR=SUMW2V/SUMW**2 VAR=MAX(VAR,WMSD) SIGF=SQRT(VAR) C C INVERT INVERSE SCALING FACTOR. C SIGF=SIGF/F**2 F=1/F C C ACCUMULATE SUMS FOR AVERAGE SCALING FACTOR AND AVERAGE C FRACTIONAL ERROR, AND TABULATE EXTREME VALUES. C NSUM=NSUM+1 SUMF=SUMF+F SUME=SUME+SIGF/F FMIN=MIN(FMIN,F) FMAX=MAX(FMAX,F) EMIN=MIN(EMIN,(SIGF/F)) EMAX=MAX(EMAX,(SIGF/F)) C C APPLY SCALING FACTOR TO Y = I/LP VALUE. C VAR=(Y*SIGF)**2+F**2*(SIGY**2+(P*Y)**2) SIGFY=SQRT(VAR) FY=F*Y C C OUTPUT C CALL WRITE1 (IO2,II,IH,IK,IL,A1,A2,A3,A4,FY,SIGFY,XTIME) NWRIT=NWRIT+1 IF (II.LT.0) GO TO 1 J=ABS(IH) K=ABS(IK) L=ABS(IL) IF (J.EQ.0.AND.K.EQ.0) GO TO 39 IF (J.EQ.0.AND.L.EQ.0) GO TO 39 IF (K.EQ.0.AND.L.EQ.0) GO TO 39 IF (J.EQ.K.AND.L.EQ.0) GO TO 39 IF (J.EQ.L.AND.K.EQ.0) GO TO 39 IF (K.EQ.L.AND.J.EQ.0) GO TO 39 IF (J.EQ.K.AND.K.EQ.L) GO TO 39 GO TO 1 39 WRITE (ILP,310) II,IH,IK,IL,Y,SIGY,FY,SIGFY,XTIME,F,SIGF NPRINT=NPRINT+1 IF (MOD(NPRINT,50).NE.0) GO TO 1 WRITE (ILP,309) ATIME,ADATE,TITLE GO TO 1 9 CONTINUE FAV=SUMF/NSUM EAV=SUME/NSUM WRITE (ILP,311) FMIN,FMAX,FAV,EMIN,EMAX,EAV IF (NWRIT.NE.NREAD) WRITE (ILP,312) NREAD,NWRIT,NREAD-NWRIT RETURN 300 FORMAT('1'/'0PROGRAM TSCALE. ',A,1X,A,'. ',A) 400 FORMAT('0LATTICE PARAMETERS:'/'0 A B C', &' ALPHA BETA GAMMA'/3F10.4,3F10.3/'0RECIPROCAL SPAC', &'E METRIC TENSOR:'/'0G11 G12 G13 ',3F10.6/' G21 G22 G23 = ', &3F10.6/' G31 G23 G33 ',3F10.6) 401 FORMAT('0FOR EACH REFLECTION, AVERAGED SCALE FACTOR WILL BE WEIGHT &ED ANISOTROPICALLY.'/'0W1 = 1/(1 + (C1*(1 - COS(PHI)**2)**2), WHER &E PHI = ANGLE(H,H0) AND'/'0C1 = ',E10.3) 402 FORMAT('0FOR EACH REFLECTION, AVERAGED SCALE FACTOR WILL BE WEIGHT &ED BY INTENSITY DIFFERENCES.'/'0W2 = 1/(1 + (C2*(I - I0)**2)'/'0C2 & = ',E10.3) 403 FORMAT('0FOR EACH REFLECTION, AVERAGED SCALE FACTOR WILL BE WEIGHT &ED BY SIN(THETA)/LAMBDA DIFFERENCES.'/'0W3 = 1/(1 + (C3*(S - S0)** &)'/'0C3 = ',E10.3) 309 FORMAT('1'/'0PROGRAM TSCALE. ',A,1X,A,'. ',A/'0 I H K ', &' L Y = I/LP SIGMA(Y) Z = K(T)*Y SIGMA(Z) XTIME = T K(T) SIGMA &(K)'/' - - - - -------- -------- ---------- -------', &'- --------- ---- --------') 310 FORMAT(' ',I6,3I4,4F10.2,F10.3,2F6.3) 319 FORMAT('+',77X,' WARNING: GAP IN ONE OR MORE SCALING CURVES.') 311 FORMAT('0MINIMUM SCALING FACTOR KMIN = ',F7.3/ & ' MAXIMUM SCALING FACTOR KMAX = ',F7.3/ & '0AVERAGE SCALING FACTOR K = ',F7.3/ & '0MINIMUM FRACTIONAL ERROR EMIN = ',F7.3/ & ' MAXIMUM FRACTIONAL ERROR EMAX = ',F7.3/ & '0AVERAGE FRACTIONAL ERROR SIGMA(K)/K = E = ',F7.3) 312 FORMAT('0NIN = ',I5,' REFLECTIONS READ FROM THE INPUT FILE.'/' NO &UT = ',I5,' REFLECTIONS WRITTEN TO THE OUTPUT FILE.'/'0NREJ = ', &I5,' REFLECTIONS REJECTED BASED ON INPUT SERIAL NUMBERS OR EXPOSUR &E TIMES,'/' OR DUE TO A GAP BETWEEN SEGMENTS OF THE SC &ALING FUNCTION.') END C----------------------------------------------------------------------- SUBROUTINE CULL C C CULLS THE STANDARD REFERENCE REFLECTIONS FROM THE INPUT C REFLECTION DATA FILE AND COMPILES SUMMARY STATISTICS ON THE C DISTRIBUTIONS OF THE REFERENCE INTENSITIES. C PARAMETER (MSTD=20,MCUT=10,MDAT=500,MOMIT=20,MREJ=50) CHARACTER TITLE*80,ATIME*8,ADATE*9 COMMON /BLOCK1/ TITLE,ATIME,ADATE COMMON /BLOCK2/ IO1,IO2,ILP,ITABL,IPLOT,IQSUM,IGRAPH, & IAPPLY,IANISO,IYDIFF,ISDIFF,A(6),GINV(3,3) COMMON /BLOCK3/ NSTD,IHKL(MSTD,3),NCUT(MSTD),XCUT(MSTD,MCUT),XZERO COMMON /BLOCK4/ NIOMIT(MSTD),IOMIT(MSTD,MOMIT,2), & NXOMIT(MSTD),XOMIT(MSTD,MOMIT,2) COMMON /BLOCK6/ NDAT(MSTD,MCUT),XX(MSTD,MCUT,MDAT), & YY(MSTD,MCUT,MDAT),SS(MSTD,MCUT,MDAT) DIMENSION SUMN(MSTD),SUMW(MSTD), & SUMY(MSTD),SUMWY(MSTD),SUMY2(MSTD),SUMWY2(MSTD), & IMIN(MSTD),TMIN(MSTD),XMIN(MSTD),YMIN(MSTD),SMIN(MSTD), & IMAX(MSTD),TMAX(MSTD),XMAX(MSTD),YMAX(MSTD),SMAX(MSTD), & H(3) DO 10 I=1,MSTD SUMN(I)=0 SUMW(I)=0 SUMY(I)=0 SUMWY(I)=0 SUMY2(I)=0 SUMWY2(I)=0 XMIN(I)=+1E6 XMAX(I)=-1E6 YMIN(I)=+1E6 YMAX(I)=-1E6 DO 10 J=1,MCUT NDAT(I,J)=0 10 CONTINUE XTIME1=+1E5 XTIME2=-1E5 IF (ITABL.EQ.0) GO TO 11 WRITE (ILP,102) ATIME,ADATE,TITLE WRITE (ILP,101) JJ=0 11 CONTINUE C C LOOP THROUGH REFLECTION DATA FILE. C REWIND IO1 IEND=0 1 CALL READ1 (IO1,IEND,II,IH,IK,IL,A1,A2,A3,A4,Y,SIGY,XTIME) IF (IEND.NE.0) GO TO 9 C C X**N IS UNDEFINED FOR X = N = 0. C IF (XTIME.EQ.0) XTIME=0.0001 XTIME1=MIN(XTIME1,XTIME) XTIME2=MAX(XTIME2,XTIME) C C IS THIS A STANDARD REFERENCE REFLECTION? C IF (II.GT.0) GO TO 1 C C LIST MEASUREMENTS OF SETS OF STANDARDS. C IF (ITABL.EQ.0) GO TO 15 IF (II.NE.JJ-1) WRITE (ILP,109) WRITE (ILP,110) II,IH,IK,IL,Y,SIGY,XTIME JJ=II 15 CONTINUE C C WHAT TYPE OF STANDARD IS IT? C IF (NSTD.EQ.0) GO TO 50 DO 20 I=1,NSTD JH=IHKL(I,1) JK=IHKL(I,2) JL=IHKL(I,3) IF (IH.NE.JH) GO TO 20 IF (IK.NE.JK) GO TO 20 IF (IL.NE.JL) GO TO 20 C C IS IT TO BE OMITTED FROM THE CALCULATION OF SCALING FACTORS? C IF (NIOMIT(I).EQ.0) GO TO 21 DO 22 J=1,NIOMIT(I) I1=IOMIT(I,J,1) I2=IOMIT(I,J,2) IF (I1.LE.ABS(II).AND.ABS(II).LE.I2) GO TO 1 22 CONTINUE 21 IF (NXOMIT(I).EQ.0) GO TO 30 DO 23 J=1,NXOMIT(I) X1=XOMIT(I,J,1) X2=XOMIT(I,J,2) IF (X1.LT.XTIME.AND.XTIME.LT.X2) GO TO 1 23 CONTINUE C C IN WHICH TIME SEGMENT IS IT? C 30 J=1 IF (NCUT(I).EQ.0) GO TO 40 DO 31 J=1,NCUT(I) IF (XTIME.LE.XCUT(I,J)) GO TO 40 31 CONTINUE J=NCUT(I)+1 GO TO 40 20 CONTINUE 50 CONTINUE C C IF IT IS A NEW TYPE, EXTEND THE TABLE OF TYPES. C NSTD=NSTD+1 I=NSTD IHKL(I,1)=IH IHKL(I,2)=IK IHKL(I,3)=IL J=1 40 CONTINUE C C TABULATE DATA FOR SUBSEQUENT ANALYSIS. C NDAT(I,J)=NDAT(I,J)+1 IF (NDAT(I,J).GT.MDAT) STOP 'TOO MANY MEASUREMENTS OF STANDARDS' K=NDAT(I,J) XX(I,J,K)=XTIME YY(I,J,K)=Y SS(I,J,K)=SIGY C C COMPILE SUMMARY STATISTICS FOR EACH TYPE OF STANDARD. C W=1/SIGY**2 SUMN(I)=SUMN(I)+1 SUMW(I)=SUMW(I)+W SUMY(I)=SUMY(I)+Y SUMWY(I)=SUMWY(I)+W*Y SUMY2(I)=SUMY2(I)+Y**2 SUMWY2(I)=SUMWY2(I)+W*Y**2 XMIN(I)=MIN(XMIN(I),XTIME) XMAX(I)=MAX(XMAX(I),XTIME) IF (YMIN(I).LE.Y) GO TO 41 IMIN(I)=II TMIN(I)=XTIME YMIN(I)=Y SMIN(I)=SIGY 41 CONTINUE IF (YMAX(I).GE.Y) GO TO 42 IMAX(I)=II TMAX(I)=XTIME YMAX(I)=Y SMAX(I)=SIGY 42 CONTINUE C C LOOP BACK FOR NEXT MEASUREMENT. C GO TO 1 9 CONTINUE C C CALCULATE AND PRINT REFERENCE TIME OF THE EXPERIMENT. C WRITE (ILP,213) XTIME1,XTIME2,0.5*(XTIME1+XTIME2) IF (XZERO.EQ.0) XZERO=0.5*(XTIME1+XTIME2) C C CHECK THAT THE DOMAIN OF EACH STANDARD INCLUDES THE REFERENCE C TIME OF THE EXPERIMENT. C DO 19 I=1,NSTD IF (SUMN(I).EQ.0) GO TO 19 IF (XMIN(I).LE.XZERO.AND.XZERO.LE.XMAX(I)) GO TO 19 WRITE (ILP,214) I,(IHKL(I,K),K=1,3),XMIN(I),XMAX(I) 19 CONTINUE C C EVALUATE AND PRINT SUMMARY STATISTICS. C WRITE (ILP,102) ATIME,ADATE,TITLE WRITE (ILP,210) DO 39 I=1,NSTD H(1)=IHKL(I,1) H(2)=IHKL(I,2) H(3)=IHKL(I,3) SINTHL=0.5*SQRT(VTMV(H,GINV,H)) N=SUMN(I) IF (N.EQ.0) GO TO 38 YMEAN=SUMY(I)/SUMN(I) WYMEAN=SUMWY(I)/SUMW(I) ESD=SQRT(SUMN(I)/SUMW(I)) IF (N.GT.1) GO TO 37 RMSD=0 RWMSD=0 GO TO 36 37 RMSD=SQRT(ABS(SUMY2(I)/SUMN(I)-YMEAN**2)*SUMN(I)/(SUMN(I)-1)) RWMSD=SQRT(ABS(SUMWY2(I)/SUMW(I)-WYMEAN**2)*SUMN(I)/(SUMN(I)-1)) 36 WRITE (ILP,211) I,(IHKL(I,K),K=1,3),SINTHL, & IMIN(I),TMIN(I),YMIN(I),SMIN(I), & N,YMEAN,RMSD,WYMEAN,RWMSD,ESD, & IMAX(I),TMAX(I),YMAX(I),SMAX(I) GO TO 39 38 WRITE (ILP,212) I,(IHKL(I,K),K=1,3),SINTHL 39 CONTINUE RETURN 102 FORMAT('1'/'0PROGRAM TSCALE. ',A,1X,A,'. ',A) 101 FORMAT ('0STANDARD REFERENCE REFLECTIONS:'/'0 J H K', &' L Y SIGY XTIME'/) 109 FORMAT (' ') 110 FORMAT (' ',I10,3X,3I3,2X,2F10.2,F10.3) 210 FORMAT('0SUMMARY STATISTICS ON STANDARD REFERENCE REFLECTIONS:'/ &'0Y = I/LP'/ &' W = 1/SIGMA(Y)**2'/ &'0MEAN = SUM(Y)/N'/ &' WMEAN = SUM(W*Y)/SUM(W)'/ &'0RMSD = SQRT((N/(N-1))*(SUM(Y**2)/N - (SUM(Y)/N)**2))'/ &' RWMSD = SQRT((N/(N-1))*(SUM(W*Y**2)/SUM(W) - (SUM(W*Y)/SUM(W))** &2))'/ &' ESD = SQRT(N/SUM(W))'/ &'0 I H K L SIN(TH)/L J T(HR) Y SIGMA(Y &) N MEAN RMSD WMEAN RWMSD ESD'/ &' ',48X,'MIN/MAX'/) 211 FORMAT(' ',I3,3X,3I3,F10.4,I10,F10.3,2F10.2,I10,5F10.2/' ',25X, &I10,F10.3,2F10.2/) 212 FORMAT(' ',I3,3X,3I3,F10.4,49X,'0'/) 213 FORMAT('0TMIN = ',F10.3,' HR'/ & ' TMAX = ',F10.3,' HR'/ & '0MID-TIME OF THE EXPERIMENT:'/ & '0TZERO = (TMIN + TMAX)/2 = ',F10.3,' HR') 214 FORMAT('0WARNING: THE DOMAIN OF THE FOLLOWING STANDARD REFLECTION & DOES NOT INCLUDE TZERO.'/' PERHAPS THE STANDARD SHOULD BE OMITTED & FROM THE CALCULATION OF SCALING FACTORS,'/' OR A DIFFERENT REFERE &NCE TIME SHOULD BE CHOSEN.'/' I H K L TMIN TMAX' &/' ',I3,3X,3I3,2F10.3) END C----------------------------------------------------------------------- FUNCTION FCRIT (NFREE) C C TABLE OF F-VALUES FOR ONE VS. NFREE DEGREES OF FREEDOM. C C PROBABILITY, P = P(F, 1, NFREE), THAT A RANDOM SAMPLE OF DATA C WOULD GIVE A RATIO OF CHI-SQUARED VALUES LARGER THAN F. C C P. R. BEVINGTON, DATA REDUCTION AND ERROR ANALYSIS FOR THE C PHYSICAL SCIENCES, MC GRAW-HILL BOOK CO., NEW YORK, 1969, C P. 318. C DIMENSION N(20),F(20) DATA N /1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, & 15, 20, 24, 30, 40, 60, 120, 1000000/ C C P = 0.1 C DATA F /39.9, 8.53, 5.54, 4.54, 4.06, 3.78, 3.59, 3.46, 3.36, & 3.28, 3.23, 3.18, 3.07, 2.97, 2.93, 2.88, 2.84, 2.79, & 2.75, 2.71/ C C P = 0.05 C C DATA F /161, 18.5, 10.1, 7.71, 6.61, 5.99, 5.59, 5.32, 5.12, C & 4.96, 4.84, 4.75, 4.54, 4.35, 4.26, 4.17, 4.08, 4.00, C & 3.92, 3.84/ C C P = 0.01 C C DATA F /4050, 98.5, 34.1, 21.2, 16.3, 13.7, 12.2, 11.3, 10.6, C & 10.0, 9.65, 9.33, 8.68, 8.10, 7.82, 7.56, 7.31, 7.08, C & 6.85, 6.63/ IF (NFREE.GT.12) GO TO 5 FCRIT=F(NFREE) RETURN 5 DO 1 I=12,19 IF (NFREE.GT.N(I+1)) GO TO 1 DF=F(I+1)-F(I) DN=N(I+1)-N(I) DFDN=DF/DN DN=NFREE-N(I) FCRIT=F(I)+DN*DFDN RETURN 1 CONTINUE END C----------------------------------------------------------------------- SUBROUTINE GRAPH C C CALCULATES AND PLOTS POINTS OF INVERSE SCALING FUNCTION C R = Y/Y0 AND ITS FRACTIONAL ERROR SIGMA(R)/R. C PARAMETER (MSTD=20,MCUT=10,MDAT=500,MOMIT=20,MREJ=50) PARAMETER (NX=101,NY=51) CHARACTER TITLE*80,ATIME*8,ADATE*9 COMMON /BLOCK1/ TITLE,ATIME,ADATE COMMON /BLOCK2/ IO1,IO2,ILP,ITABL,IPLOT,IQSUM,IGRAPH, & IAPPLY,IANISO,IYDIFF,ISDIFF,A(6),GINV(3,3) COMMON /BLOCK3/ NSTD,IHKL(MSTD,3),NCUT(MSTD),XCUT(MSTD,MCUT),XZERO COMMON /BLOCK6/ NDAT(MSTD,MCUT),XX(MSTD,MCUT,MDAT), & YY(MSTD,MCUT,MDAT),SS(MSTD,MCUT,MDAT) COMMON /BLOCK7/ PAR(MSTD,MCUT,4),COV(MSTD,MCUT,4,4), & XLIM(MSTD,MCUT,2),YZERO(MSTD),P DIMENSION X(MDAT),Y(MDAT),S(MDAT),C(4) C C FIND PLOT LIMITS. C XMIN=+1E6 XMAX=-1E6 YMIN=+1E6 YMAX=-1E6 SUMSSQ=0 N=0 DO 20 I=1,NSTD DO 20 J=1,NCUT(I)+1 IF (NDAT(I,J).EQ.0) GO TO 20 DO 21 K=1,NDAT(I,J) XI=XX(I,J,K) YI=YY(I,J,K) SIGYI=SS(I,J,K) XMIN=MIN(XMIN,XI) XMAX=MAX(XMAX,XI) YMIN=MIN(YMIN,YI-2*SIGYI) YMAX=MAX(YMAX,YI+2*SIGYI) SUMSSQ=SUMSSQ+SIGYI**2 N=N+1 21 CONTINUE 20 CONTINUE RMSSIG=SQRT(SUMSSQ/N) YMID=(YMIN+YMAX)/2 DELY=(NY-1)*RMSSIG Y1=YMID-DELY/2 Y2=YMID+DELY/2 YMIN=MIN(YMIN,Y1) YMAX=MAX(YMAX,Y2) DELX=(XMAX-XMIN)/(NX-1) DELY=(YMAX-YMIN)/(NY-1) C C PLOT DATA AND SCALING FUNCTION FOR FOR EACH STANDARD. C DO 31 I=1,NSTD N=0 DO 32 J=1,NCUT(I)+1 IF (NDAT(I,J).EQ.0) GO TO 32 DO 33 K=1,NDAT(I,J) N=N+1 X(N)=XX(I,J,K) Y(N)=YY(I,J,K) 33 CONTINUE 32 CONTINUE IF (N.EQ.0) GO TO 31 WRITE (ILP,110) ATIME,ADATE,TITLE CALL PLOT3 (I,N,X,Y,XMIN,XMAX,YMIN,YMAX) WRITE (ILP,101) (IHKL(I,K),K=1,3),DELX,DELY,RMSSIG 31 CONTINUE C C PLOT COMPOSITE SCALING FUNCTION AND ITS FRACTIONAL ERROR. C FMIN=+1E6 FMAX=-1E6 SUMF=0 DO 11 I=1,NX SUMN=0 SUMW=0 SUMWY=0 SUMWY2=0 X(I)=0 Y(I)=0 S(I)=0 XI=XMIN+(I-1)*DELX C C Q IS AN EXPANSION FACTOR USED TO DEAL WITH THE (RARE) CASE OF C A TIME X NOT IN THE DOMAIN OF ANY SEGMENT. EXTRAPOLATE NO C MORE THAN HALF THE AVERAGE TIME BETWEEN MEASUREMENTS OF C STANDARDS. C Q=0 10 IF (Q.GT.0.5) GO TO 11 DO 12 J=1,NSTD DO 12 K=1,NCUT(J)+1 IF (NDAT(J,K).EQ.0) GO TO 12 N=NDAT(J,K) X1=XLIM(J,K,1) X2=XLIM(J,K,2) DX=Q*(X2-X1)/N X1=X1-DX X2=X2+DX IF (XI.LT.X1) GO TO 12 IF (XI.GT.X2) GO TO 12 YI=0 VAR=0 DO 13 L=1,4 YI=YI+PAR(J,K,L)*XI**(L-1) DO 13 M=1,4 VAR=VAR+XI**(L+M-2)*COV(J,K,L,M) 13 CONTINUE C C VARIANCE CAN BE ZERO OR, DUE TO NUMERICAL IMPRECISION IN C COV(L,M), HAVE A SMALL NEGATIVE VALUE WHEN X IS NEAR X = XZERO C VAR=MAX(VAR,1E-9) WI=1/VAR SUMN=SUMN+1 SUMW=SUMW+WI SUMWY=SUMWY+WI*YI SUMWY2=SUMWY2+WI*YI**2 12 CONTINUE Q=Q+0.1 IF (SUMN.EQ.0) GO TO 10 C C CALCULATE MEAN INVERSE SCALING FACTOR AND ITS STANDARD C DEVIATION. C YI=SUMWY/SUMW VAR=1/SUMW WMSD=SUMWY2/SUMW-YI**2 IF (SUMN.GT.1) WMSD=WMSD/(SUMN-1) VAR=MAX(VAR,WMSD) SIG=SQRT(VAR) X(I)=XI Y(I)=YI S(I)=SIG/YI FMIN=MIN(FMIN,YI) FMAX=MAX(FMAX,YI) SUMF=SUMF+YI 11 CONTINUE FAVG=SUMF/NX C C PLOT COMPOSITE SCALING FUNCTION. C WRITE (ILP,110) ATIME,ADATE,TITLE CALL PLOT2 (NX,X,Y,XMIN,XMAX,YMIN,YMAX) WRITE (ILP,102) FAVG,FMIN,FMAX,DELX,DELY C C PLOT FRACTIONAL ERROR IN COMPOSITE SCALING FUNCTION. C SUMY=0 YMIN=0 YMAX=0 DO 14 I=1,NX YI=S(I) SUMY=SUMY+YI YMAX=MAX(YMAX,YI) 14 CONTINUE YAVG=SUMY/NX DELY=(YMAX-YMIN)/(NY-1) WRITE (ILP,110) ATIME,ADATE,TITLE CALL PLOT2 (NX,X,S,XMIN,XMAX,YMIN,YMAX) WRITE (ILP,103) YAVG,YMAX,DELX,DELY RETURN 110 FORMAT('1'/'0',10X,'PROGRAM TSCALE. ',A,1X,A,'. ',A) 101 FORMAT('0',10X,'Y = INVERSE SCALING FACTOR VERSUS X = TIME (HR).'/ & '0',10X,' H K L'/' ',10X,3I3/ & '0',10X,'DELTA(X) = ',F7.3,' HR AND DELTA(Y) = ',E10.3, & ' PER GRID POINT. R.M.S. SIGMA(Y) = 'E10.3'.') 102 FORMAT('0',10X,'Y = AVERAGE OF INVERSE SCALING FACTORS VERSUS X =' & ' TIME (HR).'/ & '0',10X,'AVERAGE SCALING FACTOR = ',F7.3,'; MINIMUM' & ' SCALING FACTOR = ',F7.3,'; MAXIMUM SCALING' & ' FACTOR = ',F7.3,'.'/ & '0',10X,'DELTA(X) = ',F7.3,' HR AND DELTA(Y) = ',E10.3, & ' PER GRID POINT.') 103 FORMAT('0',10X,'Y = FRACTIONAL ERROR IN AVERAGED INVERSE SCALING', & ' FACTOR VERSUS X = TIME (HR).'/ & '0',10X,'AVERAGE FRACTIONAL ERROR = ',F7.3,'; MAXIMUM' & ' FRACTIONAL ERROR = ',F7.3,'.'/ & '0',10X,'DELTA(X) = ',F7.3,' HR AND DELTA(Y) = ',E10.3, & ' PER GRID POINT.') END C----------------------------------------------------------------------- SUBROUTINE MATINV (N,M,A,D) C C INVERT A SQUARE MATRIX BY GAUSS-JORDAN ELIMINATION, AND CALCULATE C ITS DETERMINANT. C C A - INPUT MATRIX, WHICH IS REPLACED BY ITS INVERSE C N - ORDER OF MATRIX C - LOGICAL SIZE OF ARRAY A C M - PHYSICAL SIZE OF ARRAY A, M .GE. N C D - DETERMINANT OF THE INPUT MATRIX C C RETURNS D = 0 TO FLAG A SINGULAR MATRIX. C C MATRIX IS SCALED TO AVOID ARITHMETIC OVERFLOW OR UNDERFLOW DURING C INVERSION. C C MATRIX SCALING USES A DIAGONAL MATRIX Q WITH C C Q(I,I) = 1/SQRT(ABS(A(I,I))). C C FORM B = Q A Q, C C B(I,J) = Q(I,I)*A(I,J)*Q(J,J), C C AND INVERT B TO OBTAIN B-1. THEN, C C B-1 = (Q A Q)-1 C = (A Q)-1 Q-1 C B-1 Q = Q-1 A-1 C Q B-1 Q = A-1 C C IN ADDITION, C C DET(A) = DET(B)/PROD(I=1,N) Q(I,I)**2 C C IF N EXCEEDS MAXN, THE WORK VECTORS Q, IK, AND JK MUST BE RE- C DIMENSIONED TO A LENGTH OF N. C PARAMETER (MAXN=25) DIMENSION A(M,M),Q(MAXN),IK(MAXN),JK(MAXN) DOUBLE PRECISION A,Q,AMAX,T C C STORE RECIPROCAL SQUARE-ROOT DIAGONAL MAGNITUDES FOR SCALING. C DO 10 I=1,N T=SQRT(ABS(A(I,I))) IF (T.EQ.0) THEN D=0 RETURN ELSE Q(I)=1/T END IF 10 CONTINUE C C SCALE MATRIX. C DO 11 I=1,N DO 11 J=1,N A(I,J)=A(I,J)*Q(I)*Q(J) 11 CONTINUE C C INVERT SCALED MATRIX AND CALCULATE ITS DETERMINANT. C D=1 DO 100 K=1,N C C FIND LARGEST ELEMENT A(I,J) IN REST OF MATRIX. C AMAX=0 DO 30 I=K,N DO 30 J=K,N IF (ABS(A(I,J)).GT.ABS(AMAX)) THEN AMAX=A(I,J) IK(K)=I JK(K)=J END IF 30 CONTINUE D=D*AMAX IF (D.EQ.0) RETURN C C INTERCHANGE ROWS AND COLUMNS TO PUT AMAX IN A(K,K). C I=IK(K) IF (I.GT.K) THEN DO 50 J=1,N T=A(K,J) A(K,J)=A(I,J) A(I,J)=-T 50 CONTINUE END IF J=JK(K) IF (J.GT.K) THEN DO 60 I=1,N T=A(I,K) A(I,K)=A(I,J) A(I,J)=-T 60 CONTINUE END IF C C ACCUMULATE ELEMENTS OF INVERSE MATRIX. C DO 70 I=1,N IF (I.NE.K) A(I,K)=-A(I,K)/AMAX 70 CONTINUE DO 80 I=1,N DO 80 J=1,N IF (I.NE.K.AND.J.NE.K) A(I,J)=A(I,J)+A(I,K)*A(K,J) 80 CONTINUE DO 90 J=1,N IF (J.NE.K) A(K,J)=+A(K,J)/AMAX 90 CONTINUE A(K,K)=1/AMAX 100 CONTINUE C C RESTORE ORDERING OF MATRIX. C DO 130 K=N,1,-1 J=IK(K) IF (J.GT.K) THEN DO 110 I=1,N T=A(I,K) A(I,K)=-A(I,J) A(I,J)=T 110 CONTINUE END IF I=JK(K) IF (I.GT.K) THEN DO 120 J=1,N T=A(K,J) A(K,J)=-A(I,J) A(I,J)=T 120 CONTINUE END IF 130 CONTINUE C C SCALE INVERSE MATRIX. C DO 140 I=1,N DO 140 J=1,N A(I,J)=A(I,J)*Q(I)*Q(J) 140 CONTINUE C C SCALE DETERMINANT, IF POSSIBLE. C T=0 DO 150 I=1,N T=T+LOG(Q(I)) 150 CONTINUE T=LOG(ABS(D))-2*T C C TEST AGAINST RANGE OF MACHINE ALLOWED MAGNITUDES. C C EIGHT-BIT EXPONENT HAS MAXIMUM MAGNITUDE 2**7 = 128. C XMIN = 2**(-128) = 0.29E-38 LOG(XMIN) = -88.7 C XMAX = (2**(+128))/2 = 1.70E+38 LOG(XMAX) = +88.0 C IF (-87.LT.T.AND.T.LT.+87) D=(D/ABS(D))*EXP(T) RETURN END C----------------------------------------------------------------------- SUBROUTINE METRIC (A,GINV) C C CALCULATE RECIPROCAL METRIC TENSOR. C DIMENSION A(6),GINV(3,3) PI=ACOS(-1.0) CA=COS(A(4)*PI/180) CB=COS(A(5)*PI/180) CG=COS(A(6)*PI/180) SA=SIN(A(4)*PI/180) SB=SIN(A(5)*PI/180) SG=SIN(A(6)*PI/180) V=A(1)*A(2)*A(3)*SQRT(1-CA**2-CB**2-CG**2+2*CA*CB*CG) B1=A(2)*A(3)*SA/V B2=A(1)*A(3)*SB/V B3=A(1)*A(2)*SG/V B4=(CB*CG-CA)/(SB*SG) B5=(CA*CG-CB)/(SA*SG) B6=(CA*CB-CG)/(SA*SB) GINV(1,1)=B1*B1 GINV(1,2)=B1*B2*B6 GINV(1,3)=B1*B3*B5 GINV(2,1)=GINV(1,2) GINV(2,2)=B2*B2 GINV(2,3)=B2*B3*B4 GINV(3,1)=GINV(1,3) GINV(3,2)=GINV(2,3) GINV(3,3)=B3*B3 RETURN END C----------------------------------------------------------------------- SUBROUTINE PCALC (N,X,Y,S,C,PSQ) C C CALCULATES 'INSTABILITY CONSTANT' P AS DESCRIBED BY C MC CANDLISH, STOUT, AND ANDREWS, ACTA CRYST., A31, C 245-249 (1975). C DIMENSION X(1),Y(1),S(1),C(1) SUMDEL=0 SUMSIG=0 SUMOBS=0 DO 1 I=1,N SIGY=S(I) YOBS=Y(I) YCALC=0 DO 2 J=1,4 YCALC=YCALC+C(J)*X(I)**(J-1) 2 CONTINUE SUMDEL=SUMDEL+(YOBS-YCALC)**2 SUMSIG=SUMSIG+SIGY**2 SUMOBS=SUMOBS+YOBS**2 1 CONTINUE PSQ=(SUMDEL-SUMSIG)/SUMOBS RETURN END C----------------------------------------------------------------------- SUBROUTINE PLOTA (AA) PARAMETER (NX=101,NY=51) CHARACTER AA*1 DIMENSION AA(NX,NY) DO 10 IX=1,NX DO 10 IY=1,NY AA(IX,IY)=' ' 10 CONTINUE DO 21 IX=1,NX AA(IX, 1)='.' AA(IX,NY)='.' 21 CONTINUE DO 22 IX=1,NX,10 AA(IX, 1)='+' AA(IX,NY)='+' 22 CONTINUE DO 23 IY=1,NY AA( 1,IY)='.' AA(NX,IY)='.' 23 CONTINUE DO 24 IY=1,NY,10 AA( 1,IY)='+' AA(NX,IY)='+' 24 CONTINUE RETURN END C----------------------------------------------------------------------- SUBROUTINE PLOTB (XMIN,XMAX,YMIN,YMAX,AA) PARAMETER (NX=101,NY=51) CHARACTER AA*1 DIMENSION AA(NX,NY),XGRID(11) COMMON /BLOCK2/ IO1,IO2,ILP,ITABL,IPLOT,IQSUM,IGRAPH, & IAPPLY,IANISO,IYDIFF,ISDIFF,A(6),GINV(3,3) DELX=XMAX-XMIN DELY=YMAX-YMIN F=0 DO 90 IY=1,NY IF (MOD(IY,10).EQ.1) THEN YGRID=YMAX-F*DELY F=F+0.2 WRITE (ILP,101) YGRID,(AA(IX,IY),IX=1,NX) ELSE WRITE (ILP,100) (AA(IX,IY),IX=1,NX) END IF 90 CONTINUE 101 FORMAT(' ',E10.3,101A1) 100 FORMAT(' ',10X, 101A1) F=0 DO 91 I=1,11 XGRID(I)=XMIN+F*DELX F=F+0.1 91 CONTINUE WRITE (ILP,102) XGRID 102 FORMAT(' ',1X,11F10.2) RETURN END C----------------------------------------------------------------------- SUBROUTINE PLOT1 (N,X,Y,XMIN,XMAX,YMIN,YMAX,A) C C PLOT DATA POINTS AND FITTED CURVE FOR THE JTH SEGMENT OF THE C ITH STANDARD. C PARAMETER (NX=101,NY=51) CHARACTER AA*1 DIMENSION X(1),Y(1),A(1),AA(NX,NY) C C CLEAR THE PLOT ARRAY AND FILL IN THE GRID OUTLINE. C CALL PLOTA (AA) C C FILL IN THE CALCULATED POINTS. C DELX=XMAX-XMIN DELY=YMAX-YMIN IF (A(1).EQ.0.AND.A(2).EQ.0.AND.A(3).EQ.0.AND.A(4).EQ.0) GO TO 49 DX=DELX/(NX-1) DO 41 I=1,NX XI=XMIN+(I-1)*DX YI=0 DO 42 J=1,4 YI=YI+A(J)*XI**(J-1) 42 CONTINUE IX=1+((XI-XMIN)/DELX)*(NX-1)+0.49 IY=1+((YMAX-YI)/DELY)*(NY-1)+0.49 IF (IX.LT.1.OR.IX.GT.NX) GO TO 41 IF (IY.LT.1.OR.IY.GT.NY) GO TO 41 AA(IX,IY)='*' 41 CONTINUE 49 CONTINUE C C FILL IN THE EXPERIMENTAL DATA POINTS. C DO 30 I=1,N IX=1+((X(I)-XMIN)/DELX)*(NX-1)+0.49 IY=1+((YMAX-Y(I))/DELY)*(NY-1)+0.49 IX=MAX(IX,1) IX=MIN(IX,NX) IY=MAX(IY,1) IY=MIN(IY,NY) AA(IX,IY)='O' 30 CONTINUE C C PRINT THE PLOT ARRAY. C CALL PLOTB (XMIN,XMAX,YMIN,YMAX,AA) RETURN END C----------------------------------------------------------------------- SUBROUTINE PLOT2 (N,X,Y,XMIN,XMAX,YMIN,YMAX) C C PLOT THE Q-SUM FUNCTION, THE COMPOSITE SCALING FUNCTION, OR C THE FRACTIONAL ERROR IN THE COMPOSITE SCALING FUNCTION. C PARAMETER (NX=101,NY=51) CHARACTER AA*1 DIMENSION X(1),Y(1),AA(NX,NY) C C CLEAR THE PLOT ARRAY AND FILL IN THE GRID OUTLINE. C CALL PLOTA (AA) C C FILL IN THE CALCULATED POINTS. C DELX=XMAX-XMIN DELY=YMAX-YMIN DO 30 I=1,N IX=1+((X(I)-XMIN)/DELX)*(NX-1)+0.49 IY=1+((YMAX-Y(I))/DELY)*(NY-1)+0.49 IF (IX.LT.1.OR.IX.GT.NX) GO TO 30 IF (IY.LT.1.OR.IY.GT.NY) GO TO 30 AA(IX,IY)='*' 30 CONTINUE C C PRINT THE PLOT ARRAY. C CALL PLOTB (XMIN,XMAX,YMIN,YMAX,AA) RETURN END C----------------------------------------------------------------------- SUBROUTINE PLOT3 (ISTD,NDATA,X,Y,XMIN,XMAX,YMIN,YMAX) C C PLOT DATA POINTS AND THE (SEGMENTS OF THE) SCALING FUNCTION C FOR THE ITH STANDARD. C PARAMETER (MSTD=20,MCUT=10,MDAT=500,MOMIT=20,MREJ=50) PARAMETER (NX=101,NY=51) CHARACTER AA*1 DIMENSION X(1),Y(1),C(4),AA(NX,NY) COMMON /BLOCK2/ IO1,IO2,ILP,ITABL,IPLOT,IQSUM,IGRAPH, & IAPPLY,IANISO,IYDIFF,ISDIFF,A(6),GINV(3,3) COMMON /BLOCK3/ NSTD,IHKL(MSTD,3),NCUT(MSTD),XCUT(MSTD,MCUT),XZERO COMMON /BLOCK6/ NDAT(MSTD,MCUT),XX(MSTD,MCUT,MDAT), & YY(MSTD,MCUT,MDAT),SS(MSTD,MCUT,MDAT) COMMON /BLOCK7/ PAR(MSTD,MCUT,4),COV(MSTD,MCUT,4,4), & XLIM(MSTD,MCUT,2),YZERO(MSTD),P C C CLEAR THE PLOT ARRAY AND FILL IN THE GRID OUTLINE. C CALL PLOTA (AA) C C FILL IN THE CALCULATED POINTS. C DELX=XMAX-XMIN DELY=YMAX-YMIN DX=DELX/(NX-1) DO 40 I=1,NX XI=XMIN+(I-1)*DX C C FIND THE SEGMENT THAT INCLUDES X = XI IN ITS DOMAIN. C DO 41 JSEG=1,NCUT(ISTD)+1 X1=XLIM(ISTD,JSEG,1) X2=XLIM(ISTD,JSEG,2) IF (X1.LE.XI.AND.XI.LE.X2) GO TO 42 41 CONTINUE GO TO 40 42 CONTINUE DO 43 K=1,4 C(K)=PAR(ISTD,JSEG,K) 43 CONTINUE C C CALCULATE YI = YI(XI). C YI=0 DO 44 K=1,4 YI=YI+C(K)*XI**(K-1) 44 CONTINUE IX=1+((XI-XMIN)/DELX)*(NX-1)+0.49 IY=1+((YMAX-YI)/DELY)*(NY-1)+0.49 IF (IX.LT.1.OR.IX.GT.NX) GO TO 40 IF (IY.LT.1.OR.IY.GT.NY) GO TO 40 AA(IX,IY)='*' 40 CONTINUE C C FILL IN THE EXPERIMENTAL DATA POINTS. C DO 30 I=1,NDATA IX=1+((X(I)-XMIN)/DELX)*(NX-1)+0.49 IY=1+((YMAX-Y(I))/DELY)*(NY-1)+0.49 IX=MAX(IX,1) IX=MIN(IX,NX) IY=MAX(IY,1) IY=MIN(IY,NY) AA(IX,IY)='O' 30 CONTINUE C C PRINT THE PLOT ARRAY. C CALL PLOTB (XMIN,XMAX,YMIN,YMAX,AA) RETURN END C----------------------------------------------------------------------- SUBROUTINE INPUT C C READ CONTROL DATA. C PARAMETER (MSTD=20,MCUT=10,MDAT=500,MOMIT=20,MREJ=50) CHARACTER TITLE*80,ATIME*8,ADATE*9,FILE1*40,FILE2*40 COMMON /BLOCK1/ TITLE,ATIME,ADATE COMMON /BLOCK2/ IO1,IO2,ILP,ITABL,IPLOT,IQSUM,IGRAPH, & IAPPLY,IANISO,IYDIFF,ISDIFF,A(6),GINV(3,3) COMMON /BLOCK3/ NSTD,IHKL(MSTD,3),NCUT(MSTD),XCUT(MSTD,MCUT),XZERO COMMON /BLOCK4/ NIOMIT(MSTD),IOMIT(MSTD,MOMIT,2), & NXOMIT(MSTD),XOMIT(MSTD,MOMIT,2) COMMON /BLOCK5/ NIREJ,IREJ(MREJ,2), & NXREJ,XREJ(MREJ,2) IO1=1 IO2=2 ILP=3 OPEN (UNIT=IO1,NAME='tscale.dat',STATUS='OLD') READ (IO1,101) TITLE READ (IO1,101) FILE1 READ (IO1,101) FILE2 ITABL=1 IPLOT=1 IQSUM=0 IGRAPH=1 IAPPLY=1 IANIS0=0 IYDIFF=0 ISDIFF=0 A(1)=0 NSTD=0 DO 1000 I=1,MSTD NIOMIT(I)=0 NXOMIT(I)=0 NCUT(I)=0 1000 CONTINUE NIREJ=0 NXREJ=0 ICUT=0 IIOMIT=0 IXOMIT=0 IIREJ=0 IXREJ=0 XZERO=0 READ (IO1,104,END=9) & ITABL,IPLOT,IQSUM,IGRAPH,IAPPLY,IANISO,IYDIFF,ISDIFF C C FOR ANISOTROPIC SCALING OR SIN(THETA)/LAMBDA DEPENDENT C SCALING, READ LATTICE PARAMETERS. C READ (IO1,105,END=9) A C C TABULATE (1) ANY STANDARDS THAT ARE TO BE ANALYZED IN SEG- C MENTS, (2) AND (3) ANY STANDARDS THAT ARE TO BE OMITTED FROM C THE ANALYSIS, AND (4) AND (5) ANY DATA THAT ARE TO BE REJECTED C RATHER THAN SCALED. C C EACH OF THE FIVE INPUT LISTS IS ENDED BY A BLANK RECORD. C 70 READ (IO1,103,END=9) JH,JK,JL,XTCUT IF (JH.EQ.0.AND.JK.EQ.0.AND.JL.EQ.0) GO TO 79 ICUT=1 CALL TABULATE (JH,JK,JL,I) NCUT(I)=NCUT(I)+1 IF (NCUT(I).EQ.MCUT) STOP 'TOO MANY EXPOSURE TIME CUTS' J=NCUT(I) C C X-RAY EXPOSURE TIMES TO THE NEAREST 0.001 HR = 3.6 SEC C XCUT(I,J)=XTCUT+0.00049 GO TO 70 79 CONTINUE 10 READ (IO1,102,END=9) JH,JK,JL,I1,I2 IF (JH.EQ.0.AND.JK.EQ.0.AND.JL.EQ.0) GO TO 19 IF (I2.EQ.0) I2=I1 IIOMIT=1 CALL TABULATE (JH,JK,JL,I) NIOMIT(I)=NIOMIT(I)+1 IF (NIOMIT(I).GT.MOMIT) STOP 'TOO MANY IOMIT RANGES' J=NIOMIT(I) IOMIT(I,J,1)=ABS(I1) IOMIT(I,J,2)=ABS(I2) GO TO 10 19 CONTINUE 20 READ (IO1,103,END=9) JH,JK,JL,X1,X2 IF (JH.EQ.0.AND.JK.EQ.0.AND.JL.EQ.0) GO TO 29 IF (X2.EQ.0) X2=X1 X1=X1-0.00051 X2=X2+0.00051 IXOMIT=1 CALL TABULATE (JH,JK,JL,I) NXOMIT(I)=NXOMIT(I)+1 IF (NXOMIT(I).GT.MOMIT) STOP 'TOO MANY XOMIT RANGES' J=NXOMIT(I) XOMIT(I,J,1)=X1 XOMIT(I,J,2)=X2 GO TO 20 29 CONTINUE 30 READ (IO1,104,END=9) I1,I2 IF (I1.EQ.0.AND.I2.EQ.0) GO TO 39 IF (I2.EQ.0) I2=I1 IIREJ=1 NIREJ=NIREJ+1 IF (NIREJ.GT.MREJ) STOP 'TOO MANY IREJECT RANGES' J=NIREJ IREJ(J,1)=ABS(I1) IREJ(J,2)=ABS(I2) GO TO 30 39 CONTINUE 40 READ (IO1,105,END=9) X1,X2 IF (X1.EQ.0.AND.X2.EQ.0) GO TO 49 IF (X2.EQ.0) X2=X1 C C X-RAY EXPOSURE TIMES TO THE NEAREST 0.001 HR = 3.6 SEC C X1=X1-0.00051 X2=X2+0.00051 IXREJ=1 NXREJ=NXREJ+1 IF (NXREJ.GT.MREJ) STOP 'TOO MANY XREJECT RANGES' J=NXREJ XREJ(J,1)=X1 XREJ(J,2)=X2 GO TO 40 49 CONTINUE C C READ IN A VALUE OF XZERO, THE REFERENCE TIME FOR SCALING. C IF THE VALUE IS XZERO = 0, THE PROGRAM WILL DEFAULT TO THE C MID-TIME OF THE EXPERIMENT, XZERO = 0.5*(XMIN + XMAX). C READ (IO1,105,END=9) XZERO 9 CLOSE (UNIT=IO1) OPEN (UNIT=IO1,NAME=FILE1,STATUS='OLD',FORM='UNFORMATTED', & READONLY) IF (IAPPLY.NE.0) & OPEN (UNIT=IO2,NAME=FILE2,STATUS='NEW',FORM='UNFORMATTED') OPEN (UNIT=ILP,NAME='tscale.lp',STATUS='NEW',CARRIAGECONTROL= &'FORTRAN') C C PRINT CONTROL DATA. C CALL TIME (ATIME) CALL DATE (ADATE) WRITE (ILP,201) ATIME,ADATE,TITLE WRITE (ILP,202) FILE1 WRITE (ILP,203) FILE2 WRITE (ILP,205) & ITABL,IPLOT,IQSUM,IGRAPH,IAPPLY,IANISO,IYDIFF,ISDIFF IF (A(1).NE.0) THEN WRITE (ILP,212) A CALL METRIC (A,GINV) END IF IF (ICUT.EQ.0) GO TO 61 WRITE (ILP,220) DO 62 I=1,NSTD IF (NCUT(I).EQ.0) GO TO 62 WRITE (ILP,221) I,(IHKL(I,J),J=1,3),(XCUT(I,J),J=1,NCUT(I)) 62 CONTINUE 61 IF (IIOMIT.EQ.0.AND.IXOMIT.EQ.0) GO TO 51 WRITE (ILP,204) IF (IIOMIT.EQ.0) GO TO 52 DO 53 I=1,NSTD IF (NIOMIT(I).EQ.0) GO TO 53 DO 43 J=1,NIOMIT(I) WRITE (ILP,206) I,(IHKL(I,K),K=1,3),(IOMIT(I,J,K),K=1,2) 43 CONTINUE 53 CONTINUE 52 IF (IXOMIT.EQ.0) GO TO 51 DO 54 I=1,NSTD IF (NXOMIT(I).EQ.0) GO TO 54 DO 44 J=1,NXOMIT(I) WRITE (ILP,207) I,(IHKL(I,K),K=1,3),(XOMIT(I,J,K),K=1,2) 44 CONTINUE 54 CONTINUE 51 IF (IIREJ.EQ.0.AND.IXREJ.EQ.0) GO TO 55 WRITE (ILP,208) IF (IIREJ.EQ.0) GO TO 56 WRITE (ILP,209) ((IREJ(I,J),J=1,2),I=1,NIREJ) 56 IF (IXREJ.EQ.0) GO TO 55 WRITE (ILP,210) ((XREJ(I,J),J=1,2),I=1,NXREJ) 55 IF (XZERO.EQ.0) GO TO 59 WRITE (ILP,211) XZERO 59 RETURN 101 FORMAT(A) 102 FORMAT(3I4,2I10) 103 FORMAT(3I4,2F10.0) 104 FORMAT(8I10) 105 FORMAT(8F10.0) 201 FORMAT('1'/'0PROGRAM TSCALE. ',A,1X,A,'. ',A) 202 FORMAT('0INPUT REFLECTION DATA FILE TO BE SCALED: ',A) 203 FORMAT('0OUTPUT SCALED REFLECTION DATA FILE: ',A) 205 FORMAT('0JOB CONTROL VARIABLES:'/'0 ITABL IPLOT IQSUM IGRAP &H IAPPLY IANISO IYDIFF ISDIFF'/' ',8I8) 212 FORMAT('0LATTICE PARAMETERS:'/'0 A B C', &' ALPHA BETA GAMMA'/' ',3F10.4,3F10.3) 204 FORMAT('0STANDARD REFERENCE REFLECTION MEASUREMENTS TO BE OMITTED & FROM THE ANALYSIS:'/ &'0 I H K L FROM TO'/) 206 FORMAT(' ',I3,3X,3I3,2I10) 207 FORMAT(' ',I3,3X,3I3,2F10.3,' HR') 208 FORMAT('0REFLECTION MEASUREMENTS TO BE REJECTED FROM THE DATA SET: &'/ &'0 FROM TO'/) 209 FORMAT(' ',15X,2I10) 210 FORMAT(' ',15X,2F10.3,' HR') 220 FORMAT('0STANDARD REFERENCE REFLECTION DATA TO BE ANALYZED IN TIME & SEGMENTS:'/ &'0 I H K L XTCUT(HR)'/) 221 FORMAT(' ',I3,3X,3I3,10F10.3) 211 FORMAT('0INPUT VALUE OF XZERO:'/'0XZERO = ',F10.3,' HR'/'0WILL B', &'E USED INSTEAD OF THE MID-TIME OF THE EXPERIMENT AS THE REFERENCE & TIME FOR THE SCALING.') END C----------------------------------------------------------------------- SUBROUTINE READ1 (IFILE,IEND,II,IH,IK,IL,A1,A2,A3,A4,Y,SIGY,XTIME) IEND=0 READ (IFILE,END=9,ERR=9) II,IH,IK,IL,A1,A2,A3,A4,Y,SIGY,XTIME C C NICOLET (SYNTEX) TIME REGISTER OVERFLOWS AFTER (2**16 - 1) MIN = C 1092.25 HR. C IF (XTIME.LT.0) XTIME=XTIME+1092.25 RETURN 9 IEND=1 RETURN END C----------------------------------------------------------------------- SUBROUTINE REGRESS (NDAT,X,Y,SIG,PAR,COV) C C LEAST-SQUARES FIT OF A POWER SERIES POLYNOMIAL BY MULTIPLE C LINEAR REGRESSION. C C Y = A0 + A1*X + A2*X**2 + ... + AN*X**N C C RETURNS SCALED INVERSE MATRIX AS VARIANCE-COVARIANCE MATRIX. C DIMENSION X(1),Y(1),SIG(1) DIMENSION XMEAN(3),SUMWX(3),SUMWX2(3),SIGX(3),RXY(3),RXX(3,3) DIMENSION AXX(3,3),A(3),COVA0(3),PAR(4),COV(4,4) DOUBLE PRECISION AXX NMAX=3 DO 1 J=1,NMAX+1 PAR(J)=0 DO 1 K=1,NMAX+1 COV(J,K)=0 1 CONTINUE NMAX=MIN(NMAX,NDAT-2) C C CALCULATE WEIGHTED MEANS AND STANDARD DEVIATIONS OF YOBS AND C XOBS**J IN ORDER TO SEPARATE THE A0 PARAMETER AND SCALE THE C LEAST-SQUARES NORMAL MATRIX AND VECTOR. C SUMW=0 SUMWY=0 SUMWY2=0 DO 10 J=1,NMAX SUMWX(J)=0 SUMWX2(J)=0 10 CONTINUE DO 11 I=1,NDAT WI=1/SIG(I)**2 YI=Y(I) SUMW=SUMW+WI SUMWY=SUMWY+WI*YI SUMWY2=SUMWY2+WI*YI**2 DO 11 J=1,NMAX XIJ=X(I)**J SUMWX(J)=SUMWX(J)+WI*XIJ SUMWX2(J)=SUMWX2(J)+WI*XIJ**2 11 CONTINUE YMEAN=SUMWY/SUMW SIGY=SQRT(ABS((SUMWY2/SUMW)-YMEAN**2)*NDAT/(NDAT-1)) DO 19 J=1,NMAX XMEAN(J)=SUMWX(J)/SUMW SIGX(J)=SQRT(ABS((SUMWX2(J)/SUMW)-XMEAN(J)**2)*NDAT/(NDAT-1)) 19 CONTINUE C C ACCUMULATE NORMAL MATRIX RXX(NMAX,NMAX) AND VECTOR RXY(NMAX). C DO 30 J=1,NMAX RXY(J)=0 DO 30 K=1,J RXX(J,K)=0 30 CONTINUE DO 35 I=1,NDAT WI=1/SIG(I)**2 YI=Y(I) DO 35 J=1,NMAX XIJ=X(I)**J RXY(J)=RXY(J)+WI*(XIJ-XMEAN(J))*(YI-YMEAN) DO 35 K=1,J XIK=X(I)**K RXX(J,K)=RXX(J,K)+WI*(XIJ-XMEAN(J))*(XIK-XMEAN(K)) 35 CONTINUE DO 39 J=1,NMAX RXY(J)=(RXY(J)/SUMW)*NDAT/(NDAT-1) RXY(J)=RXY(J)/(SIGX(J)*SIGY) DO 39 K=1,J RXX(J,K)=(RXX(J,K)/SUMW)*NDAT/(NDAT-1) RXX(J,K)=RXX(J,K)/(SIGX(J)*SIGX(K)) RXX(K,J)=RXX(J,K) 39 CONTINUE C C TRY POLYNOMIALS OF SUCESSIVELY LOWER DEGREE AND PERFORM THE F- C TEST ON CHI-SQUARED VALUES FOR THE SIGNIFICANCE OF EACH HIGHER C DEGREE TERM. C CHISQ=9E9 DO 5 NPAR=NMAX,1,-1 CHISQ1=CHISQ C C INVERT THE NORMAL MATRIX RXX(NPAR,NPAR). C DO 21 J=1,NPAR DO 21 K=1,NPAR AXX(J,K)=RXX(J,K) 21 CONTINUE CALL MATINV (NPAR,NMAX,AXX,DET) IF (DET.EQ.0) RETURN C C CALCULATE PARAMETERS. C A0=YMEAN DO 40 J=1,NPAR A(J)=0 DO 41 K=1,NPAR A(J)=A(J)+RXY(K)*AXX(J,K) 41 CONTINUE A(J)=A(J)*SIGY/SIGX(J) A0=A0-A(J)*XMEAN(J) 40 CONTINUE C C CALCULATE CHI-SQUARED. C CHISQ=0 DO 50 I=1,NDAT W=1/SIG(I)**2 YOBS=Y(I) YCALC=A0 DO 51 J=1,NPAR YCALC=YCALC+A(J)*X(I)**J 51 CONTINUE CHISQ=CHISQ+W*(YOBS-YCALC)**2 50 CONTINUE C C TEST SIGNIFICANCE OF HIGHER DEGREE TERM BY F-TEST USING C TABLE FUNCTION FCRIT. C NFREE=NDAT-NPAR-1 FTEST=(CHISQ-CHISQ1)/(CHISQ1/NFREE) C C IF F-TEST PASSED, ACCEPT HIGHER DEGREE POLYNOMIAL. C IF (FTEST.GE.FCRIT(NFREE)) GO TO 6 5 CONTINUE C C TRY A CONSTANT, ZEROTH DEGREE POLYNOMIAL, Y = A0 = YMEAN. C NPAR=0 CHISQ1=CHISQ CHISQ=0 DO 20 I=1,NDAT CHISQ=CHISQ+((Y(I)-YMEAN)/SIG(I))**2 20 CONTINUE NFREE=NDAT-1 FTEST=(CHISQ-CHISQ1)/(CHISQ1/NFREE) IF (FTEST.LT.FCRIT(NFREE)) GO TO 99 6 CONTINUE C C EVALUATE FINAL POLYNOMIAL. C NPAR=NPAR+1 DO 121 J=1,NPAR DO 121 K=1,NPAR AXX(J,K)=RXX(J,K) 121 CONTINUE CALL MATINV (NPAR,NMAX,AXX,DET) IF (DET.EQ.0) RETURN C C CALCULATE PARAMETERS. C A0=YMEAN DO 140 J=1,NPAR A(J)=0 DO 141 K=1,NPAR A(J)=A(J)+RXY(K)*AXX(J,K) 141 CONTINUE A(J)=A(J)*SIGY/SIGX(J) A0=A0-A(J)*XMEAN(J) 140 CONTINUE C C CALCULATE CHI-SQUARED. C CHISQ=0 DO 150 I=1,NDAT W=1/SIG(I)**2 YOBS=Y(I) YCALC=A0 DO 151 J=1,NPAR YCALC=YCALC+A(J)*X(I)**J 151 CONTINUE CHISQ=CHISQ+W*(YOBS-YCALC)**2 150 CONTINUE C C CALCULATE PARAMETER ERROR ESTIMATES. C NFREE=NDAT-NPAR-1 CHISQ=CHISQ/NFREE DO 89 J=1,NPAR DO 89 K=1,J COV(J,K)=CHISQ*(AXX(J,K)/SUMW)*NDAT/(NDAT-1) COV(J,K)=COV(J,K)/(SIGX(J)*SIGX(K)) COV(K,J)=COV(J,K) 89 CONTINUE C C SET UP PARAMETER AND VARIANCE-COVARIANCE ARRAYS TO INCLUDE THE C A0 ELEMENTS FOR RETURN TO THE CALLING PROGRAM. C DO 90 I=(NPAR+1),2,-1 PAR(I)=A(I-1) DO 90 J=I,2,-1 COV(I,J)=COV(I-1,J-1) COV(J,I)=COV(I,J) 90 CONTINUE PAR(1)=A0 VARA0=0 DO 91 J=1,NPAR COVA0(J)=0 DO 91 K=1,NPAR COVA0(J)=COVA0(J)-XMEAN(K)*COV(K+1,J+1) VARA0=VARA0+XMEAN(J)*XMEAN(K)*COV(J+1,K+1) 91 CONTINUE COV(1,1)=VARA0 DO 92 I=1,NPAR J=1+I COV(1,J)=COVA0(I) COV(J,1)=COV(1,J) 92 CONTINUE RETURN 99 CONTINUE C C RETURN A CONSTANT, ZEROTH DEGREE POLYNOMIAL, Y = A0. C PAR(1)=YMEAN COV(1,1)=MAX(1/SUMW,((SUMWY2/SUMW)-YMEAN**2)/(NDAT-1)) RETURN END C----------------------------------------------------------------------- SUBROUTINE TABULATE (JH,JK,JL,I) PARAMETER (MSTD=20,MCUT=10,MDAT=500,MOMIT=20,MREJ=50) COMMON /BLOCK3/ NSTD,IHKL(MSTD,3),NCUT(MSTD),XCUT(MSTD,MCUT),XZERO IF (NSTD.EQ.0) GO TO 1 DO 2 I=1,NSTD IH=IHKL(I,1) IK=IHKL(I,2) IL=IHKL(I,3) IF (JH.NE.IH) GO TO 2 IF (JK.NE.IK) GO TO 2 IF (JL.NE.IL) GO TO 2 GO TO 3 2 CONTINUE 1 CONTINUE NSTD=NSTD+1 IF (NSTD.GT.MSTD) STOP 'TOO MANY TYPES OF STANDARDS' I=NSTD IHKL(I,1)=JH IHKL(I,2)=JK IHKL(I,3)=JL 3 RETURN END C----------------------------------------------------------------------- FUNCTION VTMV(U,G,V) C C BILINEAR FORM = (VECTOR TRANSPOSE)*(MATRIX)*(VECTOR) = SCALAR C BILINEAR FORM = (ROW MATRIX)*(SQUARE MATRIX)*(COLUMN MATRIX) C DIMENSION U(3),G(3,3),V(3) VTMV=0 DO 1 I=1,3 DO 1 J=1,3 VTMV=VTMV+U(J)*G(I,J)*V(I) 1 CONTINUE RETURN END C----------------------------------------------------------------------- SUBROUTINE WRITE1 (IFILE,II,IH,IK,IL,A1,A2,A3,A4,Y,SIGY,XTIME) WRITE (IFILE) II,IH,IK,IL,A1,A2,A3,A4,Y,SIGY,XTIME RETURN END