$! To split file: $ @MODEM.TXT $ CREATE PE07A.FOR C ## PE07A 19/05/70 C NAME PE07A(R) CHECK FUNCTION PE07A (M,A,B,C,X) C###### 19/05/70 LAST LIBRARY UPDATE DIMENSION A(1),B(1),C(1) IF(M-1)4,5,6 5 PE07A=(X-A(1))*C(2)+C(1) RETURN 4 PE07A=C(1) RETURN 6 J=M-1 DPC=C(M+1) SPC=C(M) DO 16 I=1,J K=M+1-I DNC=DPC*(X-A(K))+SPC SPC=-DPC*B(K)+C(K-1) DPC=DNC 16 CONTINUE PE07A=DPC*(˜À(1))+SPC RETURN END $ CREATE SPSCAL.FOR SUBROUTINE SPSCAL(X,NLOW,NHI) C C TO SCALE A FILE TO MAX. 9000. C J.C. PHILLIPS, EMBL OUTSTATION, HAMBURG. C 17/10/79 C DIMENSION X(1) FMAX=0. DO 1 I=NLOW,NHI IF(FMAX.LT.ABS(X(I))) FMAX=ABS(X(I)) 1 CONTINUE SCALE=9000./FMAX DO 2 I=NLOW,NHI X(I)=SCALE*X(I) 2 CONTINUE RETURN END $ CREATE TG01B.FOR C ## TG01B 04/08/70 C NAME TG01B(R) CHECK REAL FUNCTION TG01B*4(IX,N,U,S,D,X) C###### 04/08/70LAST LIBRARY UPDATE C C********************************************************************** C C TG01B - FUNCTION ROUTINE TO EVALUATE A CUBIC SPLINE GIVEN SPLINE C VALUES AND FIRST DERIVATIVE VALUES AT THE GIVEN KNOTS. C C THE SPLINE VALUE IS DEFINED AS ZERO OUTSIDE THE KNOT RANGE,WHICH C IS EXTENDED BY A ROUNDING ERROR FOR THE PURPOSE. C C F = TG01B(IX,N,U,S,D,X) C C IX ALLOWS CALLER TO TAKE ADVANTAGE OF SPLINE PARAMETERS SET C ON A PREVIOUS CALL IN CASES WHEN X POINT FOLLOWS PREVIOUS C X POINT. IF IX < 0 THE WHOLE RANGE IS SEARCHED FOR KNOT C INTERVAL; IF IX > 0 IT IS ASSUMED THAT X IS GREATER THAN C THE X OF THE PREVIOUS CALL AND SEARCH STARTED FROM THERE. C N THE NUMBER OF KNOTS. C U THE KNOTS. C S THE SPLINE VALUES. C D THE FIRST DERIVATIVE VALUES OF THE SPLINE AT THE KNOTS. C X THE POINT AT WHICH THE SPLINE VALUE IS REQUIRED. C F THE VALUE OF THE SPLINE AT THE POINT X. C C MODIFIED JULY 1970 C C********************************************************************** C C ALLOWABLE ROUNDING ERROR ON POINTS AT EXTREAMS OF KNOT RANGE C IS 2**IEPS*MAX(!U(1)!,!U(N)!). INTEGER*4 IFLG,IEPS DIMENSION U(1),S(1),D(1) DATA IFLG/0/,IEPS/-19/ C C TEST WETHER POINT IN RANGE. IF(X.LT.U(1)) GO TO 990 IF(X.GT.U(N)) GO TO 991 C C JUMP IF KNOT INTERVAL REQUIRES RANDOM SEARCH. IF(IX.LT.0.OR.IFLG.EQ.0) GO TO 12 C JUMP IF KNOT INTERVAL SAME AS LAST TIME. IF(X.LE.U(J+1)) GO TO 8 C LOOP TILL INTERVAL FOUND. 1 J=J+1 11 IF(X.GT.U(J+1)) GO TO 1 GO TO 7 C C ESTIMATE KNOT INTERVAL BY ASSUMING EQUALLY SPACED KNOTS. 12 J=ABS(X-U(1))/(U(N)-U(1))*(N-1)+1 C ENSURE CASE X=U(N) GIVES J=N-1. J=MIN0(J,N-1) C INDICATE THAT KNOT INTERVAL INSIDE RANGE HAS BEEN USED. IFLG=1 C SEARCH FOR KNOT INTERVAL CONTAINING X. IF(X.GE.U(J)) GO TO 11 2 J=J-1 IF(X.LT.U(J)) GO TO 2 C C CALCULATE SPLINE PARAMETERS FOR JTH INTERVAL. 7 H=U(J+1)-U(J) Q1=H*D(J) Q2=H*D(J+1) SS=S(J+1)-S(J) B=3.0*SS-2.0*Q1-Q2 A=Q1+Q2-2.0*SS C C CALCULATE SPLINE VALUE. 8 Z=(X-U(J))/H TG01B=((A*Z+B)*Z+Q1)*Z+S(J) RETURN C TEST IF X WITHIN ROUNDING ERROR OF U(1). 990 IF(X.LE.U(1)-2.0**IEPS*AMAX1(ABS(U(1)),ABS(U(N)))) GO TO 99 J=1 GO TO 7 C TEST IF X WITHIN ROUNDING ERROR OF U(N). 991 IF(X.GE.U(N)+2.0**IEPS*AMAX1(ABS(U(1)),ABS(U(N)))) GO TO 99 J=N-1 GO TO 7 99 IFLG=0 C FUNCTION VALUE SET TO ZERO FOR POINTS OUTSIDE THE RANGE. TG01B=0.0 RETURN END $ CREATE VB06A.FOR SUBROUTINE VB06A (M,N,XD,YD,WD,RD,XN,FN,GN,DN,THETA,IPRINT,W) LOGICAL SWITCH DIMENSION XD(1),YD(1),WD(1),RD(1),XN(1),FN(1),GN(1),DN(1), 1THETA(1),W(1) C RESERVE THE FIRST FIFTEEN LOCATIONS OF W FOR C COMPUTED VALUES AND THIRD DERIVATIVES OF B-SPLINES W(1)=0. W(7)=0. C THEN SET THE KNOT POSITIONS INCLUDING SIX OUTSIDE THE RANGE C AND DELETE ANY KNOTS NOT IN STRICTLY ASCENDING ORDER I=1 NN=1 10 W(NN+18)=XN(I) 20 IF (I.GE.N) GO TO 30 I=I+1 IF (XN(I).LE.W(NN+18)) GO TO 20 NN=NN+1 XN(NN)=XN(I) THETA(NN)=THETA(I) GO TO 10 30 IF (NN.GE.N) GO TO 50 N=NN WRITE(5,40) 40 FORMAT (//5X,'SOME KNOTS HAVE BEEN DELETED BY VB06A') 50 IF (NN.GE.2) GO TO 70 WRITE(5,60) 60 FORMAT (//5X,'THE EXECUTION OF VB06A HAS BEEN STOPPED ', 1'BECAUSE THERE ARE TOO FEW KNOTS') GO TO 770 70 DO 80 I=1,3 W(19-I)=2.*W(20-I)-W(21-I) 80 W(NN+I+18)=2.*W(NN+I+17)-W(NN+I+16) C DELETE ANY DATA POINTS THAT ARE NOT IN ASCENDING ORDER I=1 MM=1 XD(MM)=AMIN1(AMAX1(XD(1),XN(1)),XN(NN)) 90 IF (I.GE.M) GO TO 100 I=I+1 IF (XD(I).LT.XD(MM)) GO TO 90 MM=MM+1 XD(MM)=AMIN1(XD(I),XN(NN)) YD(MM)=YD(I) WD(MM)=WD(I) GO TO 90 100 IF (MM.GE.M) GO TO 120 M=MM WRITE(5,110) 110 FORMAT (//5X,'SOME DATA POINTS HAVE BEEN DELETED BY VB06A') C COMPUTE THE THIRD DERIVATIVES OF THE B-SPLINES AT XN(1) 120 SWITCH=.TRUE. IN=1 GO TO 280 130 DO 140 I=2,5 W(I+6)=W(I) 140 W(I+10)=W(I) C INITIALIZE THE FACTORIZATION OF THE LEAST SQUARES MATRIX KS=NN+22 KU=KS+22 DO 150 I=KS,KU 150 W(I)=0. IM=1 IN=2 KK=KS SWITCH=.FALSE. C TEST WHETHER A DATA POINT OR A KNOT DEFINES THE NEXT EQUATION 160 KR=KK NC=5 IF (IM.GT.MM) GO TO 260 IF (XD(IM).GT.XN(IN)) GO TO 260 C CALCULATE THE VALUES OF THE B-SPLINES AT XD(IM) 170 W(2)=1./(W(IN+18)-W(IN+17)) DO 190 J=2,4 I=J+1 W(I)=0. 180 II=IN+I+16 W(I)=((XD(IM)-W(II-J))*W(I-1)+(W(II)-XD(IM))*W(I))/(W(II)-W(II-J)) I=I-1 IF (I.GE.2) GO TO 180 190 CONTINUE IF (SWITCH) GO TO 540 C COMPLETE THE EQUATION OF THE DATA POINT W(6)=YD(IM) C=WD(IM)**2 IM=IM+1 C REVISE NC PRIOR TO THE NEXT GIVENS TRANSFORMATION 200 NC=NC-1 IF (NC.LE.0) GO TO 160 C MAKE THE GIVENS TRANSFORMATION 210 ALPHA=C*W(2) IF (ALPHA.NE.0.) GO TO 230 DO 220 I=1,NC 220 W(I+1)=W(I+2) GO TO 250 230 DSTAR=W(KR)+ALPHA*W(2) BETA=W(KR)/DSTAR C=BETA*C GAMMA=ALPHA/DSTAR W(KR)=DSTAR DSTAR=W(2) DO 240 I=1,NC W(I+1)=W(I+2)-DSTAR*W(KR+I) 240 W(KR+I)=BETA*W(KR+I)+GAMMA*W(I+2) 250 KR=KR+7 GO TO 200 C ADVANCE THE CURRENT ROW OF THE UPPER TRIANGULAR MATRIX 260 IF (IN.GE.NN) GO TO 330 W(KK+28)=0. DO 270 I=4,28,6 W(KK+I+1)=W(KK+I) 270 W(KK+I)=0. KK=KK+7 C CALCULATE THE THIRD DERIVATIVES OF THE B-SPLINES AT XN(IN+) 280 W(2)=6./(W(IN+19)-W(IN+18)) DO 300 J=2,4 I=J+1 W(I)=0. 290 II=IN+I+17 W(I)=(W(I-1)-W(I))/(W(II)-W(II-J)) I=I-1 IF (I.GE.2) GO TO 290 300 CONTINUE IF (SWITCH) GO TO 130 C SET UP THE EQUATION FOR THE SMOOTHING TERM AT A KNOT ALPHA=0. I=6 310 W(I)=W(I-1)-ALPHA I=I-1 IF (I.LE.1) GO TO 320 ALPHA=W(I+6) W(I+6)=W(I) GO TO 310 320 C=THETA(IN)**2 IN=IN+1 GO TO 210 C BACK-SUBSTITUTE TO OBTAIN THE MULTIPLIERS OF THE B-SPLINES 330 KU=KK KUU=KU+21 KB=KS+7 IBS=5 DO 335 I=3,6 KK=KK+7 335 W(KK-2)=W(KK-I) 340 KK=KUU+IBS KR=KK 350 KK=KK-7 K=KK J=KK-IBS 360 J=J+1 K=K+7 W(KK)=W(KK)-W(J)*W(K) IF (K.LT.KR) GO TO 360 KR=MIN0(KR,KK+21) IF (KK.GT.KB) GO TO 350 IF (IBS-6) 370,470,500 370 SWITCH=.TRUE. ALPHA=1. DO 380 K=KS,KUU,7 380 ALPHA=AMIN1(ALPHA,W(K)) IF (ALPHA.LE.0.) GO TO 592 IF (NN.LE.2) GO TO 600 C NEXT A SPLINE WITH FIXED VALUES OF S'''(XN(1)) AND S'''(XN(N)) C IS CALCULATED TO OBTAIN THE OPTIMAL EXTREME VALUES OF S'''(X). C SET THE COEFFICIENTS OF S'''(X(1)) AND S'''(X(N)) KR=KS KK=KU DO 410 I=8,11 W(KK+4)=W(I) KK=KK+7 W(KR+6)=W(I+4) 410 KR=KR+7 DO 420 K=KR,KUU,7 420 W(K+6)=0. IBS=6 KRR=KS C BACK-SUBSTITUTE USING THE TRANSPOSE OF THE UPPER TRIANGULAR C MATRIX TO CALCULATE THE RESIDUALS OF EACH EXTRA SPLINE 425 I=KRR+7 KR=KRR DO 450 KK=I,KUU,7 J=KK K=KK 430 K=K-7 J=J-6 W(KK+IBS)=W(KK+IBS)-W(J)*W(K+IBS) IF (K.GT.KR) GO TO 430 450 KR=MAX0(KR,KK-21) DO 460 KK=KRR,KUU,7 460 W(KK+IBS)=W(KK+IBS)/W(KK) IF (IBS.EQ.6) GO TO 340 DO 465 K=KS,KUU,7 W(K+7)=0. 465 IF (K.GE.KU) W(K+7)=W(K+4) IBS=7 GO TO 340 C BACK-SUBSTITUTE TO GIVE THE VECTOR FOR S'''(X(N)) 470 IBS=4 KRR=KU GO TO 425 C AT EACH DATA POINT CALCULATE THE RESIDUALS OF THREE SPLINES 500 DO 520 I=6,10 520 W(I)=0. IN=2 KK=KS DO 580 IM=1,MM 530 IF (XD(IM).LE.XN(IN)) GO TO 170 IN=IN+1 KK=KK+7 GO TO 530 540 KR=KK W(11)=YD(IM) W(12)=0. W(13)=0. DO 560 I=2,5 DO 550 J=5,7 550 W(J+6)=W(J+6)-W(I)*W(KR+J) 560 KR=KR+7 DO 570 I=11,13 570 W(I)=WD(IM)*W(I) C FORM THE NORMAL EQUATIONS TO FIX S'''(X) AT THE ENDS OF THE RANGE DO 580 I=12,13 K=I+I-29 DO 580 J=11,I 580 W(K+J)=W(K+J)+W(I)*W(J) C CALCULATE THE B-SPLINE MULTIPLIERS OF THE REQUIRED FIT ALPHA=(W(8)*W(9)-W(6)*W(10))/(W(7)*W(10)-W(9)**2) BETA=(-W(8)-ALPHA*W(9))/W(10) DO 590 K=KS,KUU,7 590 W(K+5)=W(K+5)+ALPHA*W(K+6)+BETA*W(K+7) GO TO 600 C PRINT A DIAGNOSTIC MESSAGE IF THERE IS INSUFFICIENT DATA 592 WRITE(5,594) 594 FORMAT (//5X,'THE RESULTS FROM VB06A ARE UNRELIABLE BECAUSE ', 1'THERE IS INSUFFICIENT DATA TO DEFINE THE SPLINE UNIQUELY') C CALCULATE THE REQUIRED VALUES OF S(X) AND S'(X) AT THE KNOTS 600 KK=KS DO 620 IN=1,NN FN(IN)=0. GN(IN)=0. W(2)=(W(IN+19)-W(IN+18))/((W(IN+19)-W(IN+17))*(W(IN+19)-W(IN+16))) W(3)=(W(IN+18)-W(IN+17))/((W(IN+19)-W(IN+17))*(W(IN+20)-W(IN+17))) W(4)=0. DO 610 I=1,3 J=IN+I+14 ALPHA=((XN(IN)-W(J))*W(I)+(W(J+4)-XN(IN))*W(I+1))/(W(J+4)-W(J)) FN(IN)=FN(IN)+ALPHA*W(KK+5) BETA=3.*(W(I)-W(I+1))/(W(J+4)-W(J)) GN(IN)=GN(IN)+BETA*W(KK+5) 610 KK=KK+7 620 KK=KK-14 C CALCULATE THE THIRD DERIVATIVE DISCONTINUITIES OF THE SPLINE DN(1)=0. IN=1 IM=1 630 H=XN(IN+1)-XN(IN) ALPHA=FN(IN)-FN(IN+1)+H*GN(IN) BETA=ALPHA+FN(IN)-FN(IN+1)+H*GN(IN+1) DN(IN+1)=6.*BETA/H**3 DN(IN)=DN(IN+1)-DN(IN) IN=IN+1 C CALCULATE THE RESIDUALS AT THE DATA POINTS 640 IF (IM.GT.MM) GO TO 650 IF (XD(IM).GT.XN(IN)) GO TO 630 C=(XD(IM)-XN(IN-1))/H RD(IM)=YD(IM)-C*FN(IN)+(C-1.)*(FN(IN-1)+C*(ALPHA-C*BETA)) IM=IM+1 GO TO 640 650 IF (IN.LT.NN) GO TO 630 C PROVIDE PRINTING IF REQUESTED IF (IPRINT.LE.0) GO TO 770 WRITE(5,660) 660 FORMAT (1H1,35X,'SPLINE APPROXIMATION OBTAINED BY VB06A'//4X,'I', 19X,'XN(I)',18X,'FN(I)',18X,'GN(I)',13X,'3RD DERIV CHANGE', 211X,'THETA(I)'//) I=1 670 WRITE(5,680) I,XN(I),FN(I),GN(I) 680 FORMAT (I5,5E23.14) 690 I=I+1 IF (I-NN) 700,670,710 700 WRITE(5,680) I,XN(I),FN(I),GN(I),DN(I),THETA(I) GO TO 690 710 WRITE(5,720) 720 FORMAT (///4X,'I',9X,'XD(I)',18X,'YD(I)',18X,'WD(I)',19X,'FIT', 118X,'RESIDUAL'//) IN=2 DO 760 IM=1,MM 730 IF (XD(IM).LE.XN(IN)) GO TO 750 IN=IN+1 IF (SWITCH) GO TO 730 SWITCH=.TRUE. WRITE(5,740) 740 FORMAT (5X) GO TO 730 750 C=YD(IM)-RD(IM) WRITE(5,680) IM,XD(IM),YD(IM),WD(IM),C,RD(IM) 760 SWITCH=.FALSE. 770 RETURN END $ CREATE VC01A.FOR CC## VC01A 25/10/74 CC NAME VC01A(R) CHECK SUBROUTINE VC01A(X,Y,W,Z,N,A,B,C,G,H,L,M) C###### 25/10/74 LAST LIBRARY UPDATE COMMON U DIMENSION X(1),Y(1),W(1),A(1),B(1),C(1),G(1),H(1),L(1) DIMENSION U(1),Z(1),E(20),F(20) 4 FORMAT(3X,I5,5E16.6,6X,I5) LP=0 EJ=0.0 FJ=0.0 GJ=0.0 DO 1 I=1,N IV=I+N U(I)=SQRT(W(I)) Z(I)=Y(I)*U(I) EJ=EJ+X(I)*U(I)**2 FJ=FJ+U(I)*Z(I) GJ=GJ+U(I)**2 U(IV)=0.0 1 CONTINUE A(1)=EJ/GJ B(1)=0.0 C(1)=FJ/GJ G(1)=1.0/GJ DO 2 J=1,M EJ=0.0 FJ=0.0 GJ=0.0 H1=0.0 L2=0 DO 3 I=1,N Z(I)= Z(I)-C(J) *U(I) IV=I+N WS=(X(I)-A(J))*U(I)-B(J)*U(IV) U(IV)=U(I) U(I)=WS EJ=EJ+X(I)*U(I)**2 FJ=FJ+U(I)*Z(I) GJ=GJ+U(I)**2 H1=H1+Z(I)**2 IF(1-I)13,3,3 13 L2=(SIGN(1.0,-Z(I)*Z(I-1))+1.0)/2.0+L2 3 CONTINUE A(J+1)=EJ/GJ B(J+1)=GJ*G(J) C(J+1)=FJ/GJ G(J+1)=1.0/GJ H(J)=H1 L(J)=L2 2 CONTINUE 5 M1=M+1 IF(LP.GT.0)WRITE(LP,20) 20 FORMAT(1H1,8X,4HX(I),12X,4HY(I),12X,4HZ(I),//) H1=0.0 L2=0 DO 10 I=1,N W2=0.0 YY=C(M+1) DO 101 J=1,M J1=M-J+1 W1=W2 W2=YY YY=C(J1)+(X(I)-A(J1))*W2-B(J1+1)*W1 101 CONTINUE Z(I)=YY U(I)=(Z(I)-Y(I)) H1=H1+U(I)*U(I)*W(I) IF(1-I)111,110,110 111 L2=L2+(SIGN(1.0,-U(I)*U(I-1))+1.0)/2.0 110 IF(LP.GT.0)WRITE(LP,11)X(I),Y(I),Z(I) 10 CONTINUE H(M1)=H1 L(M1)=L2 H1=H1/(N-M1) DO 21 J=1,M1 G(J)=G(J)*H1 21 CONTINUE IF(LP.LE.0)GO TO 40 IF(LP.GT.0)WRITE(LP,7) 7 FORMAT(1H1,61X,8HVARIANCE,10X,9HRESIDUALS) IF(LP.GT.0)WRITE(LP,6) 6 FORMAT(7X,1HJ,7X,4HA(J),12X,4HB(J),12X,4HC(J),11X,7HOF C(J),6X, 114HSUM OF SQUARES,3X,12HSIGN CHANGES) M1=M+1 DO 24 I=1,M1 II=I-1 24 WRITE(LP,4)II,A(I),B(I),C(I),G(I),H(I),L(I) 11 FORMAT(3E16.6) 40 RETURN END BLOCK DATA COMMON/VC01B/LP INTEGER LP DATA LP /6/ END $ CREATE VTPLOT.FOR SUBROUTINE VTPLOT(IX,IY,MODE) C C J.C. PHILLIPS SUNY/BNL 2/10/83 C C TO UTILISE RETROGRAPHICS CAPABILITIES OF THE VT100 C TO PLOT INCOMING DATA. C BYTE BUFFER(10),BOX(16) INTEGER IPARAM(6) DATA BOX /"40,"140,"40,"100,"64,"176,"40,"100, Q"64,"176,"77,"137,"40,"140,"77,"137/ GOTO(10,20,20,20,20,10) MODE 10 CALL GETADR(IPARAM(1),BUFFER(1)) ! DRAW BOX BUFFER(1)="33 BUFFER(2)="14 BUFFER(3)=' ' BUFFER(4)="35 IF(MODE.EQ.6) BUFFER(4)="30 IPARAM(2)=4 CALL QIO("410,5,1,,,IPARAM) IF(MODE.EQ.6) RETURN DO 1 I=1,13,4 CALL GETADR(IPARAM(1),BOX(I)) CALL QIO("410,5,1,,,IPARAM) 1 CONTINUE CALL GETADR(IPARAM(1),BOX(1)) CALL QIO("410,5,1,,,IPARAM) 20 BUFFER(4)=(IX/"40).OR."40 BUFFER(5)=(IX.AND."37).OR."100 BUFFER(2)=(IY/"40).OR."40 BUFFER(3)=(IY.AND."37).OR."140 IF(MODE.NE.2) CALL GETADR(IPARAM(1),BUFFER(1)) IF(MODE.EQ.2) CALL GETADR(IPARAM(1),BUFFER(2)) IPARAM(2)=6 BUFFER(1)="35 BUFFER(6)="37 IF(MODE.EQ.2) IPARAM(2)=4 IF(MODE.EQ.3) IPARAM(2)=5 CALL QIO("410,5,1,,,IPARAM) RETURN END $ CREATE WFILE.FOR SUBROUTINE WFILE(X,FNAM1,NPNTS) DIMENSION X(512) BYTE CHAR,FNAM1(10) DATA CHAR /' '/ DATA JUNK /0/ CALL ASSIGN(1,FNAM1,10) WRITE(1,102) CHAR WRITE(1,102) CHAR WRITE(1,103) NPNTS WRITE(1,103) JUNK WRITE(1,104) (X(I), I=1,NPNTS) CALL CLOSE(1) RETURN 102 FORMAT(Q,A1) 103 FORMAT(3X,I5) 104 FORMAT(8F10.2) END $ CREATE RFILE.FOR SUBROUTINE RFILE(X,FNAM1,NPNTS) DIMENSION X(512) BYTE CHAR,FNAM1(10) CALL ASSIGN(1,FNAM1,10) READ(1,102) CHAR READ(1,102) CHAR READ(1,103) NPNTS READ(1,103) JUNK READ(1,104) (X(I), I=1,NPNTS) CALL CLOSE(1) RETURN 102 FORMAT(Q,A1) 103 FORMAT(3X,I5) 104 FORMAT(8F10.2) END $ CREATE XYPLOT.FOR SUBROUTINE XYPLOT(X,Y,NLOW,NHI,IPEN,XMIN,XPER,YMIN,YPER,ICROSS) C C J.C. PHILLIPS 10/24/83 C C TO UTILISE RETROGRAPHICS CAPABILITIES OF THE VT100 C DIMENSION X(512),Y(512),XCROSS(10),YCROSS(10) INTEGER NPCR(10) BYTE BUFFER(6),BOX(16),SCROSS(10) INTEGER IPARAM(6) COMMON /HAIRS/ XCROSS,YCROSS,NCROSS,NPCR,SCROSS DATA BOX /"40,"140,"40,"100,"64,"176,"40,"100, Q"64,"176,"77,"137,"40,"140,"77,"137/ NCROSS=0 CALL GETADR(IPARAM(1),BUFFER(1)) NPNTS=NHI-NLOW+1 IF(ABS(XPER).GT.0.00001) GO TO 2 WRITE(5,100) ! CLEAR SCREEN ETC. BUFFER(1)="33 BUFFER(2)="14 BUFFER(3)="7 BUFFER(4)="35 IPARAM(2)=4 CALL QIO("410,5,1,,,IPARAM) XPER=1.1*(X(NHI)-X(NLOW)) ! CALCULATE MAXIMUM AND MINIMUM XMIN=X(NLOW)-XPER/22. YMIN=Y(NLOW) YPER=YMIN DO 20 I=NLOW,NHI IF(YMIN.GT.Y(I)) YMIN=Y(I) IF(YPER.LT.Y(I)) YPER=Y(I) 20 CONTINUE C DRAW BOX DO 1 I=1,13,4 CALL GETADR(IPARAM(1),BOX(I)) CALL QIO("410,5,1,,,IPARAM) 1 CONTINUE CALL GETADR(IPARAM(1),BOX(1)) CALL QIO("410,5,1,,,IPARAM) CALL GETADR(IPARAM(1),BUFFER(1)) IPARAM(2)=6 BUFFER(1)="35 BUFFER(2)="64 BUFFER(3)="176 BUFFER(4)="40 BUFFER(5)="100 BUFFER(6)="37 CALL QIO("410,5,1,,,IPARAM) WRITE(5,102) NHI,X(NHI),YPER,NLOW,X(NLOW),YMIN YPER=1.1*(YPER-YMIN) YMIN=YMIN-YPER/22. 2 BUFFER(1)="35 BUFFER(2)="33 IPARAM(2)=1 IF(XPER.LT.0.0001.OR.YPER.LT.0.0001) GO TO 900 GOTO(11,12,13,14,15,16) IPEN ! SET PEN 12 BUFFER(1)="34 GO TO 11 13 BUFFER(3)="141 GO TO 17 14 BUFFER(3)="142 GO TO 17 15 BUFFER(3)="143 GO TO 17 16 BUFFER(3)="144 17 IPARAM(2)=3 11 CALL QIO("410,5,1,,,IPARAM) IPARAM(2)=4 DO 3 I=NLOW,NHI ! PLOT DATA IX=IFIX((X(I)-XMIN)*1023./XPER) BUFFER(3)=(IX/"40).OR."40 BUFFER(4)=(IX.AND."37).OR."100 IY=IFIX((Y(I)-YMIN)*670./YPER) BUFFER(1)=(IY/"40).OR."40 BUFFER(2)=(IY.AND."37).OR."140 CALL QIO("410,5,1,,,IPARAM) 3 CONTINUE NCROSS=0 4 IF(ICROSS.NE.1.OR.NCROSS.EQ.10) GO TO 5 NCROSS=NCROSS+1 BUFFER(1)="35 BUFFER(2)="33 BUFFER(3)="32 IPARAM(2)=3 CALL QIO("410,5,1,,,IPARAM) READ(5,101) SCROSS(NCROSS),(BUFFER(I), I=1,4) IX=(BUFFER(1).AND."37) IY=(BUFFER(2).AND."37) IX=32*IX+IY XCROSS(NCROSS)=XMIN+XPER*IX/1023. DO 6 I=NLOW,NHI IF(X(I).LT.XCROSS(NCROSS)) GO TO 6 NPCR(NCROSS)=I GO TO 8 6 CONTINUE 8 IY=(BUFFER(3).AND."37) IX=(BUFFER(4).AND."37) IY=32*IY+IX YCROSS(NCROSS)=YMIN+YPER*IY/670. IF(SCROSS(NCROSS).NE."121) GO TO 4 NCROSS=NCROSS-1 900 CONTINUE 5 BUFFER(1)="35 BUFFER(2)="70 BUFFER(3)="153 BUFFER(4)="40 BUFFER(5)="100 BUFFER(6)="37 IPARAM(2)=6 CALL QIO("410,5,1,,,IPARAM) IF(XPER.GT.0.0001.AND.YPER.GT.0.0001) RETURN WRITE(5,103) RETURN 100 FORMAT(30(/)) 101 FORMAT(5A1) 102 FORMAT(' ',40X,I3,2(1X,E15.8),28(/),' ',5X,I3,2(1X,E15.8)) 103 FORMAT(' SCALES TOO SMALL NO DATA PLOTTED') END $ CREATE DIR104.DIR 16-SEP-85 19:10 EXAFS.FTN;27 2. 07-DEC-84 14:00 SLICE.FTN;42 8. 07-DEC-84 14:00 VTPLOT.FTN;5 3. 24-JAN-85 18:27 REDUCE.FTN;46 11. 14-SEP-85 14:10 COMM2.FTN;10 3. 07-DEC-84 14:00 DIGIT.FTN;106 6. 07-DEC-84 14:00 PREFIT.FTN;12 10. 30-AUG-85 21:26 EXINIT.FTN;101 5. 22-MAR-85 14:27 SCPLOT.FTN;23 4. 07-DEC-84 14:01 SCTEST.FTN;14 1. 07-DEC-84 14:01 EDSPEC.FTN;3 5. 10-MAR-85 15:11 EXCOMM.FTN;126 18. 07-DEC-84 14:01 TEKPLT.FTN;1 14. 07-DEC-84 14:01 RFILE.FTN;2 1. 07-DEC-84 14:01 STEPS.FTN;10 3. 12-JUN-85 22:07 RCOMP.FTN;1 3. 07-DEC-84 14:01 FILSRT.FTN;1 2. 07-DEC-84 14:01 REPCAL.FTN;1 4. 07-DEC-84 14:01 READET.FTN;22 2. 07-DEC-84 14:01 BACOFF.FTN;4 19. 11-JAN-85 17:12 EXMOVE.FTN;36 5. 14-MAY-85 23:02 SPEAK.FTN;2 1. 04-FEB-85 14:47 GAPCAL.FTN;2 1. 07-DEC-84 14:01 CAMAC.FTN;17 5. 07-DEC-84 14:01 LTFNA.FTN;1 2. 07-DEC-84 14:01 COLEC.FTN;13 8. 07-DEC-84 14:01 PLOTS.FTN;1 1. 07-DEC-84 14:01 CINDY.FTN;5 12. 07-DEC-84 14:01 BRMOVE.FTN;7 3. 07-DEC-84 14:01 NEWFIL.FTN;1 1. 07-DEC-84 14:01 RSMOVE.FTN;11 2. 07-DEC-84 14:01 ALPHA.FTN;1 1. 07-DEC-84 14:01 CAIOX.FTN;4 4. 07-DEC-84 14:01 PRCALC.FTN;6 13. 12-JUN-85 22:14 PLOT.FTN;23 5. 25-JAN-85 16:02 WFILE.FTN;2 1. 07-DEC-84 14:02 REFINE.FTN;17 11. 11-SEP-85 10:32 SPSUM.FTN;4 4. 07-DEC-84 14:02 BAKOFF.FTN;12 20. 03-SEP-85 16:21 EXSCAN.FTN;154 14. 14-MAY-85 22:48 REDUC2.FTN;6 11. 29-APR-85 14:43 XYPLOT.FTN;13 7. 25-JAN-85 17:26 VTTEST.FTN;2 1. 07-DEC-84 14:02 TALKTEST.FTN;11 2. 07-DEC-84 14:02 REDUC3.FTN;3 11. 29-APR-85 14:47 COMMON.FTN;13 2. 28-JAN-85 09:46 SPAVER.FTN;22 10. 29-AUG-85 19:18 LINEAR.FTN;3 2. 07-DEC-84 14:02 EXPLOT.FTN;57 8. 28-JAN-85 10:07 EXREAD.FTN;35 2. 07-DEC-84 14:02 DATYPE.FTN;10 2. 07-DEC-84 14:02 REDUC4.FTN;1 11. 14-SEP-85 14:12 VERIFY.FTN;4 1. 07-DEC-84 14:03 TEST2.FTN;1 1. 07-DEC-84 14:03 UNLOCK.FTN;6 1. 07-DEC-84 14:03 SCAN.FTN;10 3. 07-DEC-84 14:03 FASCAN.FTN;51 7. 07-DEC-84 14:03 COUNT2.FTN;22 3. 07-DEC-84 14:04 FILMAK.FTN;27 16. 07-DEC-84 14:04 COUNT.FTN;46 6. 07-DEC-84 14:04 EXNORM.FTN;3 3. 07-DEC-84 14:04 DATCON.FTN;1 2. 07-DEC-84 14:04 DRAW.FTN;7 3. 07-DEC-84 14:04 REDRAW.FTN;10 2. 07-DEC-84 14:04 QCOUNT.FTN;13 3. 07-DEC-8Lš:04 DITEST.FTN;4 2. 07-DEC-84 14:05 ANOUNC.FTN;4 2. 07-DEC-84 14:05 REALEC.FTN;12 3. 07-DEC-84 14:05 FMAC3.FTN;3 1. 07-DEC-84 14:05 CEE.FTN;1 1. 07-DEC-84 14:05 CTEST3.FTN;3 2. 07-DEC-84 14:05 RELIN.FTN;7 1. 07-DEC-84 14:05 SEMOVE.FTN;34 7. 07-DEC-84 14:05 VTCOPL.FTN;4 3. 02-JAN-85 15:34 XYCOPL.FTN;21 8. 02-JAN-85 15:41 IN4027.FTN;3 1. 02-JAN-85 14:53 BUFILL.FTN;2 1. 02-JAN-85 14:44 VECTOR.FTN;22 1. 02-JAN-85 14:46 MOBEAM.FTN;2 1. 02-JAN-85 14:55 TOTAL OF 392./417. BLOCKS IN 79. FILES $ CREATE DIR103.DIR 16-SEP-85 19:11 FDR2.FTN;1 11. 07-DEC-84 11:52 PHASE.FTN;1 6. 07-DEC-84 11:52 PPOLD.FTN;1 4. 07-DEC-84 11:52 BELL.FTN;1 1. 07-DEC-84 11:52 PP.FTN;1 13. 07-DEC-84 11:52 PPDATA.FTN;1 13. 07-DEC-84 11:53 FPP2.FTN;1 2. 07-DEC-84 11:53 PE07A.FTN;1 4. 07-DEC-84 11:53 PE08A.FTN;1 7. 07-DEC-84 11:53 FORMAT.FTN;1 2. 07-DEC-84 11:53 VICTOR.FTN;1 5. 07-DEC-84 11:53 PCALC.FTN;1 2. 07-DEC-84 11:53 GRAPH.FTN;1 4. 07-DEC-84 11:53 MASS.FTN;1 3. 07-DEC-84 11:53 BSPLIN1.FTN;1 10. 07-DEC-84 11:53 MULMUL.FTN;1 3. 07-DEC-84 11:53 EXAFDR.FTN;1 11. 07-DEC-84 11:53 VC01A.FTN;1 5. 07-DEC-84 11:53 FDRSP.FTN;1 12. 07-DEC-84 11:53 BACK.FTN;1 3. 07-DEC-84 11:53 FPPFIT.FTN;1 2. 07-DEC-84 11:53 MAX.FTN;1 4. 07-DEC-84 11:53 VC03A.FTN;1 39. 07-DEC-84 11:53 CONVRT.FTN;1 3. 07-DEC-84 11:53 CCALC.FTN;1 1. 07-DEC-84 11:53 FRACT.FTN;1 3. 07-DEC-84 11:53 BACK2.FTN;1 7. 07-DEC-84 11:53 ADCAR.FTN;1 12. 07-DEC-84 11:53 AIRY.FTN;1 3. 07-DEC-84 11:53 BSUB2.FTN;1 8. 07-DEC-84 11:53 MAVER.FTN;1 4. 07-DEC-84 11:53 SCALE.FTN;1 4. 07-DEC-84 11:53 BACK3.FTN;1 7. 07-DEC-84 11:53 BSPLIN.FTN;1 11. 07-DEC-84 11:53 ANGCAL.FTN;1 2. 07-DEC-84 11:53 TG01B.FTN;1 13. 07-DEC-84 11:53 SCANEW.FTN;1 3. 07-DEC-84 11:53 EXDIAG.FTN;1 9. 07-DEC-84 11:53 SPSUM.FTN;1 4. 07-DEC-84 11:53 ESHIFT.FTN;1 3. 07-DEC-84 11:53 VB06A.FTN;1 44. 07-DEC-84 11:53 KK.FTN;1 9. 07-DEC-84 11:53 SYNANG.FTN;1 12. 07-DEC-84 11:53 ASAS.FTN;1 3. 07-DEC-84 11:53 MPLOT.FTN;1 5. 07-DEC-84 11:53 CLOSE.FTN;1 1. 07-DEC-84 11:53 PRTEST.FTN;1 11. 07-DEC-84 11:53 SPSCAL.FTN;1 1. 07-DEC-84 11:53 EXNORM.FTN;1 4. 07-DEC-84 11:53 PARFIL.FTN;1 2. 07-DEC-84 11:53 XSPLIN.FTN;1 4. 07-DEC-84 11:53 FOURT.FTN;1 4. 07-DEC-84 11:53 FIREAD.FTN;1 1. 07-DEC-84 11:53 ADRED.FTN;1 4. 07-DEC-84 11:53 FORMER.FTN;1 3. 07-DEC-84 11:53 EXFLDR.FTN;1 11. 07-DEC-84 11:53 DATEX.FTN;1 23. 07-DEC-84 11:53 EXAPOL.FTN;1 5. 07-DEC-84 11:53 PHASIN.FTN;1 2. 07-DEC-84 11:53 SPDER.FTN;1 3. 07-DEC-84 11:53 EDSPEC.FTN;1 5. 07-DEC-84 11:53 INTPOL.FTN;1 4. 07-DEC-84 11:53 STREAD.FTN;1 3. 07-DEC-84 11:53 KNIGHT.FTN;1 4. 07-DEC-84 11:53 STEPS.FTN;1 4. 07-DEC-84 11:53 SHORT.FTN;1 2. 07-DEC-84 11:53 SUBTST.FTN;1 1. 07-DEC-84 11:53 MOSTEP.FTN;1 3. 07-DEC-84 11:53 MUCALC.FTN;1 4. 07-DEC-84 11:53 AVSPEC.FTN;1 4. 07-DEC-84 11:53 BSUB.FTN;1 6. 07-DEC-84 11:53 MSPEC.FTN;2 1. 18-JUN-85 18:02 KNIN.FTN;3 1. 07-DEC-84 11:53 LIFTIM.FTN;2 2. 22-FEB-85 14:22 HUBRED.FTN;21 5. 18-JUL-85 13:07 RAINPL.FTN;12 3. 07-DEC-84 11:54 KKTRAN.FTN;2 8. 07-DEC-84 11:54 FPPCAL.FTN;10 4. 07-DEC-84 11:54 FTOTAL.FTN;2 2. 07-DEC-84 11:54 TOTAL OF 481./485. BLOCKS IN 79. FILES $ CREATE FOURT.FOR C PROGRAM FOURT C C TO FOURIER TRANSFORM A SPECTRUM. C C J.C. PHILLIPS, EMBL OUTSTATION, HAMBURG. C 9/12/79 C DIMENSION FK(1024),F(1024),FP(1024),R(256),FI(256),X(256),Y(256) INTEGER*4 INDIC(20),ISPEC BYTE HEAD(80),FILNAM(10),FNAM2(10) ISPEC=-1 IFN=10 TWOPI=6.2318 C C READ IN INPUT FILENAMES, NO. OF POINTS,RANGE. C WRITE(5,100) READ(5,101)FILNAM CALL RDSTF(1,FILNAM,IFN,ISPEC,HEAD,HEAD,INDIC,FK,NCHAN) WRITE(5,102) READ(5,101) FILNAM CALL RDSTF(1,FILNAM,IFN,ISPEC,HEAD,HEAD,INDIC,F,NCHAN) INDIC(1)=256 10 WRITE(5,105) READ(5,106) NLOW,NHI,NK,RANGE IF(NLOW.LE.0) NLOW=1 IF(NHI.LE.0) CALL EXIT IF(RANGE.LT.0.1) RANGE=5.0 DELX=RANGE/256. FKMID=(FK(NLOW)-FK(NHI))/2. SIGSQ=(FKMID**2)/2. C C ENTER OUTPUT FILENAMES. C WRITE(5,103) READ(5,101) FILNAM WRITE(5,104) READ(5,101) FNAM2 NTOP=NHI-1 NBOT=NLOW+1 DO 1 I=NBOT,NTOP FP(I)=F(I)*(FK(I+1)-FK(I-1))*EXP(-((FK(I)-FKMID)**2)/SIGSQ) IF(NK.GT.0) FP(I)=FP(I)*(FK(I)**NK) 1 CONTINUE CLOW=F(NLOW)*(FK(NBOT)-FK(NLOW)) CLOW=CLOW*EXP(-((FK(NLOW)-FKMID)**2)/SIGSQ) CHI=F(NHI)*(FK(NHI)-FK(NTOP))*EXP(-((FK(NHI)-FKMID)**2)/SIGSQ) IF(NK.LT.1) GO TO 2 CLOW=CLOW*(FK(NLOW)**NK) CHI=CHI*(FK(NHI)**NK) 2 DO 3 I=1,256 X(I)=DELX*I*TWOPI R(I)=CLOW*COS(FK(NLOW)*X(I))+CHI*COS(FK(NHI)*X(I)) FI(I)=CLOW*SIN(FK(NLOW)*X(I))+CHI*SIN(FK(NHI)*X(I)) DO 4 J=NBOT,NTOP R(I)=R(I)+FP(J)*COS(FK(J)*X(I)) FI(I)=FI(I)+FP(J)*SIN(FK(J)*X(I)) 4 CONTINUE X(I)=X(I)/TWOPI Y(I)=SQRT(R(I)**2+FI(I)**2) 3 CONTINUE C C OUTPUT RESULTS C CALL SPSCAL(Y,1,256) CALL WRSTF(2,FILNAM,IFN,ISPEC,HEAD,HEAD,INDIC,X,256) CALL WRSTF(2,FNAM2,IFN,ISPEC,HEAD,HEAD,INDIC,Y,256) GO TO 10 100 FORMAT('$ K FILENAME ? ') 101 FORMAT(10A1) 102 FORMAT('$ F(K) FILENAME ? ') 103 FORMAT('$ X FILENAME ? ') 104 FORMAT('$ F(X) FILENAME ? ') 105 FORMAT('$ NLOW,NHI,NK,RANGE ? ') 106 FORMAT(3I5,F12.5) END PIP>TT1:=EXNORM.FTN C PROGRAM EXNORM C C TO NORMALISE EXAFS SPECTRA FOR FALLOFF OF ATOMIC ABSORPTION C C J.C. PHILLIPS, EMBL OUTSTATION, HAMBURG C MODIFIED 2/22/83 FOR USE AT BNL. C DIMENSION X(512),Y(512) BYTE FILNAM(10),FNAM2(10) C C READ IN CONSTANTS AND FILENAMES C WRITE(5,111) READ(5,112) WL,C1,D1,A1 WRITE(5,100) READ(5,103) FILNAM CALL RFILE(X,FILNAM,NCH1) WRITE(5,102) READ(5,103) FNAM2 WRITE(5,104) READ(5,105) NSPEC IF(NSPEC.EQ.0) NSPEC=1 WRITE(5,107) READ(5,103) FILNAM EDGE=C1*WL**3-D1*WL**4+A1 EZERO=12398./WL DO 10 I=1,NCH1 X(I)=.0381002*X(I)**2+EZERO X(I)=12398./X(I) X(I)=C1*X(I)**3-D1*X(I)**4+A1 X(I)=EDGE/X(I) 10 CONTINUE DO 20 J=1,NSPEC CALL RFILE(Y,FNAM2,NCH2) WRITE(5,109) FNAM2 DO 3 I=1,NCH2 Y(I)=Y(I)*X(I) 3 CONTINUE CALL WFILE(Y,FILNAM,NCH2) WRITE(5,110) FILNAM CALL NEWFIL(FILNAM) CALL NEWFIL(FNAM2) 20 CONTINUE CALL EXIT 100 FORMAT('$ K-SPACE FILENAME ? ') 102 FORMAT('$ 1ST INPUT FILENAME ? ') 103 FORMAT(10A1) 104 FORMAT('$ NO. OF SPECTRA ? ') 105 FORMAT(I5) 107 FORMAT('$ 1ST OUTPUT FILENAME ? ') 110 FORMAT(' ',10A1,/) 109 FORMAT(' ',10A1) 111 FORMAT('$ EDGE WL, C1, D1, A1 ? ') 112 FORMAT(4F12.5) END $ CREATE GAGLAS.PAR 20 0 -35000.00 3000.00 3500.00 15000.00 16000.00 1000.00 5471.00 -945.70 23.75 31370.00 4000.00 2050.00 0.05 1079.00 409.60 1898.00 1675.00 -354.20 6.75 0.00