File: DINGLE.CD of Tape: Various/ETH/cd1
(Source file text)
DINGLE,9621,CM60000,CT15.
PERMF,LGO.
FTN(BL,OPT,R).
CATALOG,LGO,DINGBN.
.EOR.
PROGRAM DINGLE(OUTPUT=240B,TAPE4=OUTPUT,TAPE5,TAPE6)
C
COMMON /IN/ ANF(6,6),AF(6,6),ETA(6,6),DAF(6,6),DETA(6,6)
$ ,DK,HANG,HI,NP,T(6),DT(6),MR,ADENT(6),MF,FREQ(6)
COMMON /MOD/ AM(6,6),A0(6,6),A0M(6),DA0M(6),AF0(6,6)
$ ,ANFM(6),ETA0(6,6)
COMMON /LINE/ DY(6),A,B,DA,DB,SA,SB,SAB,RR,RMSD,CHISQ
C
DIMENSION X(400),Y(400),NB(400)
DIMENSION XL(10),YL(10)
C
100 FORMAT(I3)
110 FORMAT((F8.1))
130 FORMAT((A6))
C
READ(5,100) MF
READ(5,110) (FREQ(I),I=1,MF)
READ(5,100) MR
READ(5,130) (ADENT(I),I=1,MR)
CALL INPUT
CALL AZERO
CALL PRINT(1,0)
N=MR
DO 50 I=1,MF
DO 51 J=1,MR
YL(J)=ETA(I,J)
DY(J)=DETA(I,J)
XL(J)=T(J)
51 CONTINUE
CALL LINREG(XL,YL,N)
CALL PRINT(2,I)
TMAX=0.
DO 10 K=1,MR
10 TMAX=AMAX1(TMAX,T(K))
LMIN=1
LMAX=(400-4*MR)/2
DO 61 L=1,LMAX
X(L)=1.+(TMAX-1.)*(L-LMIN)/FLOAT(LMAX-LMIN)
Y(L)=A*X(L)+B
NB(L)=4
61 CONTINUE
DO 52 J=1,MR
YL(J)=ETA0(I,J)
DY(J)=DETA(I,J)
XL(J)=T(J)
52 CONTINUE
CALL LINREG(XL,YL,N)
CALL PRINT(3,I)
LMIN=LMAX+1
LMAX=LMAX*2
DO 62 L=LMIN,LMAX
X(L)=1.+(TMAX-1.)*(L-LMIN)/FLOAT(LMAX-LMIN)
Y(L)=A*X(L)+B
NB(L)=10
62 CONTINUE
DO 60 J=1,MR
K=4*(J-1)+LMAX+1
X(K)=T(J)
Y(K)=ETA(I,J)+DETA(I,J)
NB(K)=2
X(K+1)=T(J)
Y(K+1)=ETA(I,J)-DETA(I,J)
NB(K+1)=2
X(K+2)=T(J)
Y(K+2)=ETA(I,J)
NB(K+2)=1
X(K+3)=T(J)
Y(K+3)=ETA0(I,J)
NB(K+3)=3
60 CONTINUE
NPOINT=400
CALL GRAPH(X,Y,NB,NPOINT)
50 CONTINUE
END
SUBROUTINE INPUT
C
COMMON /IN/ ANF(6,6),AF(6,6),ETA(6,6),DAF(6,6),DETA(6,6)
$ ,DK,HANG,HI,NP,T(6),DT(6),MR,ADENT(6),MF,FREQ(6)
C
610 FORMAT(A2)
620 FORMAT(A6)
630 FORMAT(I2)
640 FORMAT(8(G10.3))
C
ICHAN=6
DO 10 I=1,MR
1 REWIND ICHAN
2 READ(ICHAN,610) AMIST
IF (AMIST.NE.2H$$) GO TO 2
READ(ICHAN,630) IMIST
IF (IMIST.NE.1) GO TO 2
READ(ICHAN,620) ANAME
IF (ANAME.NE.ADENT(I)) GO TO 2
READ(ICHAN,630) IMIST
READ(ICHAN,640) DK,HANG,HI,AMIST,T(I),DT(I),AMIST,FNP
NP=IFIX(FNP)
DK=DK*1.E-3
READ(ICHAN,610) AMIST
READ(ICHAN,630) IMIST
READ(ICHAN,630) NFREQ
DO 20 J=1,NFREQ
READ(ICHAN,640) Q1,Q2,Q3,Q4,DQ1,DQ2,DQ3,DQ4
DO 30 K=1,MF
IF (ABS(Q1-FREQ(K)).GT.1./(DK*NP)) GO TO 30
ANF(K,I)=Q1
AF(K,I)=Q3
ETA(K,I)=Q4
DAF(K,I)=DQ3
DETA(K,I)=DQ4
30 CONTINUE
20 CONTINUE
10 CONTINUE
RETURN
END
SUBROUTINE PRINT(IFUN,LOOP)
C
COMMON /IN/ ANF(6,6),AF(6,6),ETA(6,6),DAF(6,6),DETA(6,6)
$ ,DK,HANG,HI,NP,T(6),DT(6),MR,ADENT(6),MF,FREQ(6)
COMMON /MOD/ AM(6,6),A0(6,6),A0M(6),DA0M(6),AF0(6,6)
$ ,ANFM(6),ETA0(6,6)
COMMON /LINE/ DY(6),A,B,DA,DB,SA,SB,SAB,RR,RMSD,CHISQ
C
500 FORMAT(1H1)
505 FORMAT(1X,T50,#VERSION #,I3)
510 FORMAT(1X,T10,#RUNS:#,T20,A6,T30,#TO#,T40,A6,//)
520 FORMAT(1X,T20,#ANGLE= #,F6.1,T40,#HI= #,F7.3,# T#,T60
$ ,#DK= #,1PE10.2,# 1/T#,//)
530 FORMAT(1X,T20,#FREQUENCY : #,F8.3,# T#,//)
540 FORMAT(1X,T10,#RUN#,T22,#T[K]#,T32,#FREQ[T]#,T43,
$ #AMP(AF)#,T52,#VAR.AMP.#,T62,#ETA[T]#,T71,#VAR.ETA.#,
$ T80,#/#,T91,#ETA0[T]#,T105,#A0#,//)
550 FORMAT(1X,T10,A6,T20,0PF6.3,T30,F9.2,T40,1PE10.2,T50,E10.2,
$ T60,0PF7.3,T70,F7.3,T80,#/#,T90,F7.3,T100,1PE10.2)
560 FORMAT(1X,//,T10,#MEAN FREQUENCY= #,T30,F9.4,# T#,T90,
$ #MEAN A0= #,T100,1PE10.2,/,T10,#PERIOD= #,T30,1PE9.2,
$ # 1/T#,T90,#VAR.A0= #,T100,1PE10.2,////)
570 FORMAT(1X,T10,#ETA0 FIT:#,T20,#X =#,T30,F7.3,T40,#K#,T70,
$ #MU =#,T80,F7.3,/,T10,#********#,/,T33,#+/-#,T38,F7.3,
$ T83,#+/-#,T88,F7.3,# (LINE FIT)#,/,T38,F7.3,T88,F7.3,
$ # (VARIANCE)#,,//)
580 FORMAT(1X,T20,#MU*P =#,T30,3PF7.3,T40,#1E-3*1/T#,T70,
$ #CHISQ=#,T80,0PF7.3,/,T70,#RMSD=#,T80,F7.3,/,T70,
$ #RR=#,T80,F7.3,//)
590 FORMAT(1X,T10,#ETA FIT:#,T20,#X =#,T30,F7.3,T40,#K#,T70,
$ #MU =#,T80,F7.3,/,T10,#........#,/,T33,#+/-#,T38,F7.3,
$ T83,#+/-#,T88,F7.3,# (LINE FIT)#,/,T38,F7.3,T88,F7.3,
$ # (VARIANCE)#,,//)
C
C !!!!!!!!!!!!!!!!!!!!!!!
IVERS=4
C !!!!!!!!!!!!!!!!!!!!!!!
I=LOOP
IF (IFUN.EQ.2) GO TO 2
IF (IFUN.EQ.3) GO TO 3
WRITE(4,500)
WRITE(4,505) IVERS
WRITE(4,510) ADENT(1),ADENT(MR)
WRITE(4,520) HANG,HI,DK
RETURN
2 WRITE(4,500)
WRITE(4,530) FREQ(I)
WRITE(4,540)
WRITE(4,550) (ADENT(J),T(J),ANF(I,J),AF(I,J),DAF(I,J)
$ ,ETA(I,J),DETA(I,J),ETA0(I,J),A0(I,J),J=1,MR)
PER=1./ANFM(I)
WRITE(4,560) ANFM(I),A0M(I),PER,DA0M(I)
3 TD=B/A
AMU=A/14.69
DTDL=ABS(TD)*SQRT((DA/A)**2+(DB/B)**2)
TRY=SA/A**2+SB/B**2-2*SAB/A/B
IF (TRY.LE.0.) TRY=1000000.
DTDS=ABS(TD)*SQRT(TRY)
DMUL=DA/14.69
DMUS=SQRT(SA)/14.69
IF (IFUN.EQ.2) WRITE(4,590) TD,AMU,DTDL,DMUL,DTDS,DMUS
IF (IFUN.EQ.3) WRITE(4,570) TD,AMU,DTDL,DMUL,DTDS,DMUS
AMUP=AMU*PER
WRITE(4,580) AMUP,CHISQ,RMSD,RR
IF (IFUN.EQ.3) WRITE(4,500)
RETURN
END
SUBROUTINE AZERO
C
COMMON /IN/ ANF(6,6),AF(6,6),ETA(6,6),DAF(6,6),DETA(6,6)
$ ,DK,HANG,HI,NP,T(6),DT(6),MR,ADENT(6),MF,FREQ(6)
COMMON /MOD/ AM(6,6),A0(6,6),A0M(6),DA0M(6),AF0(6,6)
$ ,ANFM(6),ETA0(6,6)
C
DO 10 J=1,MR
DO 11 I=1,MF
TAU=ETA(I,J)*DK*NP
AM(I,J)=AF(I,J)*TAU/(1-EXP(-TAU))
11 A0(I,J)=AM(I,J)*SQRT(HI)*EXP(ETA(I,J)/HI)/T(J)
10 CONTINUE
DO 12 I=1,MF
SP=0.
SA=0.
SB=0.
SC=0.
DO 13 J=1,MR
A0P=A0(I,J)**2*((DAF(I,J)/AF(I,J))**2+(DETA(I,J)/ETA(I,J))**2)
A0P=1./A0P
SA=SA+A0(I,J)*A0P
SB=SB+A0(I,J)**2*A0P
SC=SC+ANF(I,J)
SP=SP+A0P
13 CONTINUE
A0M(I)=SA/SP
DA0M(I)=SQRT((SB/SP-A0M(I)**2)/(MR-1))
ANFM(I)=SC/MR
12 CONTINUE
DO 14 I=1,MF
DO 15 J=1,MR
ETA0(I,J)=HI*ALOG(A0M(I)*T(J)/SQRT(HI)/AM(I,J))
TAU0=ETA0(I,J)*DK*NP
AF0(I,J)=AM(I,J)*(1.-EXP(-TAU0))/TAU0
15 CONTINUE
14 CONTINUE
RETURN
END
SUBROUTINE LINREG(X,Y,N)
C
COMMON /LINE/ DY(6),A,B,DA,DB,SA,SB,SAB,RR,RMSD,CHISQ
C
DIMENSION V(6),P(6)
DIMENSION X(1),Y(1)
C
C LINE FIT: Y=A*X+B TO N POINTS X(I),Y(I)
C DA,DB VARIANCES FROM LINE FIT
C SA,SB VARIANCES, SAB CORRELATION FROM DY(I)
C RR=REGRESSION COEFFICIENT, RMSD=ROOT MEAN SQUARE DEVIATION
SX=0.
SY=0.
SXX=0.
SXY=0.
SYY=0.
SP=0.
DO 60 J=1,N
P(J)=1.
IF (DY(J).EQ.0.) GO TO 72
P(J)=1./DY(J)**2
72 SX=SX+X(J)*P(J)
SY=SY+Y(J)*P(J)
SXX=SXX+X(J)*X(J)*P(J)
SXY=SXY+X(J)*Y(J)*P(J)
SYY=SYY+Y(J)*Y(J)*P(J)
SP=SP+P(J)
60 CONTINUE
D=SP*SXX-SX*SX
A=(SP*SXY-SX*SY)/D
B=(SY*SXX-SX*SXY)/D
C
C
C VARIANCES FROM LINE FIT
IF (N.GT.2) GO TO 75
DA=0.
DB=0.
RR=0.
RMSD=0.
CHISQ=0.
GO TO 76
75 CONTINUE
Q=0.
QA=0.
DO 62 J=1,N
V(J)=Y(J)-X(J)*A-B
Q=Q+V(J)*V(J)*P(J)
QA=QA+(V(J))**2
62 CONTINUE
QM=(Q/(N-2))
C=SP*SYY-SY*SY
DA=SQRT(QM*SP/D)
DB=SQRT(QM*SXX/D)
RR=SQRT(A*A*D/C)
RMSD=SQRT(QA/N)
CHISQ=QM
76 CONTINUE
C VARIANCES FROM DY(I)
IF (DY(1).NE.0.) GO TO 77
SA=0.
SB=0.
SAB=0.
RETURN
77 SA=SP/D
SB=SXX/D
SAB=SX/D
RETURN
END
SUBROUTINE GRAPH(X,Y,NB,NPOINT)
C
DIMENSION X(1),Y(1),NB(1),LINEA(120),LRUB(120),XSCALE(24)
LOGICAL SQUARE
C
DATA LINEA /120*1H /
C
1000 FORMAT(3X,A6,2H :/10X,1H:)
1001 FORMAT(10X,1H:,120A1)
1002 FORMAT(1X,F8.3,2H :,120A1)
1003 FORMAT(10X,1H-,24A5)
1004 FORMAT(8X,11(F6.3,4X))
C
IXAUT=2
XMIN=1.
XMAX=-1E100
IYAUT=3
YMIN=1.E100
YMAX=-1.E100
KCOL=91.
LINE=55.
SQUARE=.F.
XAXIS=5H T[K]
YAXIS=6HETA[T]
C
LRUB(1)=1
IF(IXAUT.EQ.0) GO TO 20
DO 10 N=1,NPOINT
IF (IXAUT.NE.2) XMIN=AMIN1(XMIN,X(N))
IF (IXAUT.NE.1) XMAX=AMAX1(XMAX,X(N))
10 CONTINUE
20 IF(IYAUT.EQ.0) GO TO 40
DO 30 N=1,NPOINT
IF (IYAUT.NE.2) YMIN=AMIN1(YMIN,Y(N))
30 IF (IYAUT.NE.1) YMAX=AMAX1(YMAX,Y(N))
40 DX=(XMAX-XMIN)/(KCOL-1)
DY=(YMAX-YMIN)/(LINE-1)
IF(.NOT.SQUARE) GOTO 50
DX=DX*10.
DY=DY*6.
DX=AMAX1(DX,DY)
DY=DX
DX=DX/10.
DY=DY/6.
50 CENX=(XMIN+XMAX)/2.
CENY=(YMIN+YMAX)/2.
RKCEN=FLOAT(KCOL)/2.
RLCEN=FLOAT(LINE)/2.
WRITE(4,1000) YAXIS
NP=NPOINT
DO 100 L=1,LINE
KLAST=0
LPOINT=0
IF(NP.EQ.0) GOTO 65
DO 60 N=1,NP
NLP=N-LPOINT
X(NLP)=X(N)
Y(NLP)=Y(N)
NB(NLP)=NB(N)
IF((RLCEN-(Y(N)-CENY)/DY).GT.FLOAT(L)) GOTO 60
LPOINT=LPOINT+1
K=IFIX(RKCEN+(X(N)-CENX)/DX)+1
LRUB(LPOINT)=K
KLAST=MAX0(K,KLAST)
IF(NB(N).NE.1) GOTO 51
LINEA(K)=1HE
GOTO 60
51 IF(NB(N).NE.2) GOTO 52
LINEA(K)=1H-
GOTO 60
52 IF(NB(N).NE.3) GOTO 53
LINEA(K)=1H0
GOTO 60
53 IF(NB(N).NE.4) GOTO 54
LINEA(K)=1H.
GOTO 60
54 IF(NB(N).NE.5) GOTO 55
LINEA(K)=1H5
GOTO 60
55 IF(NB(N).NE.6) GOTO 56
LINEA(K)=1H6
GOTO 60
56 LINEA(K)=1H*
60 CONTINUE
65 LINEL=LINE-L
LSCALE=(LINEL)/6
IF((6*LSCALE).EQ.(LINEL)) GOTO 70
WRITE(4,1001) (LINEA(K),K=1,KLAST)
GOTO 80
70 YSCALE=CENY+DY*(FLOAT(LINEL)-RLCEN+.5)
WRITE(4,1002) YSCALE,(LINEA(K),K=1,KLAST)
80 NP=NP-LPOINT
DO 90 K=1,LPOINT
KL=LRUB(K)
90 LINEA(KL)=1H
100 CONTINUE
KSCALE=KCOL/5
IF(KSCALE.EQ.24) KSCALE=23
DO 110 K=1,KSCALE
110 XSCALE(K)=5H^----
WRITE(4,1003) (XSCALE(K),K=1,KSCALE),XAXIS
KSCALE=KCOL/10
IF(KSCALE.EQ.12) KSCALE=11
DO 120 K=1,KSCALE
120 XSCALE(K)=CENX+DX*(10.*FLOAT(K-1)-RKCEN+.5)
WRITE(4,1004) (XSCALE(K),K=1,KSCALE)
RETURN
END