File: DINGLE.FT of Tape: Various/ETH/prog1
(Source file text)
CDC DINGLE,9621,CM60000,CT15.
CDC PERMF,LGO.
CDC FTN(BL,R).
CDC CATALOG,LGO,DINGBN.
CDC PROGRAM DINGLE(OUTPUT=240B,TAPE4=OUTPUT,TAPE5,TAPE6)
C
DIMENSION X(400),Y(400),NB(400)
REAL MU,KET
LOGICAL ERRMU,ERR
C
100 FORMAT(I3)
110 FORMAT(F10.1)
130 FORMAT((A6))
200 FORMAT(///T8'MU NEGATIV !!!!!!!!!!!!!!'/1H1)
410 FORMAT(1X,'J= ',1I3,' X=',F7.3,' +/-',F7.3,/)
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)
READ(5,100) I0
CALL INPUT
CALL PRINT(1)
C
C
DO 10 I=1,MF
IF=I
C
DO 50 J=1,MR
AF0(J)=AF(I,J)
ETA0(J)=ETA(I,J)
TAUD=ETA(I,J)*DK*NP
DTAUD=DETA(I,J)*DK*NP
KET=(1.-EXP(-TAUD))
YL(J)=H0AL*ALOG(ABS(AF(I,J)*TAUD/T(J)/KET/I0DH))-
$I0*TAUD/NP
DY(J)=H0AL*SQRT((DAF(I,J)/AF(I,J))**2+
$((1./TAUD+EXP(-TAUD)/KET-I0/NP)*DTAUD)**2)
TDI(J)=ETA(I,J)/14.69/AMU-T(J)
TDP(J)=(DETA(I,J)/14.69/AMU)**2+(DT(J))**2
DTD=SQRT(TDP(J))
TDP(J)=1./TDP(J)
WRITE(4,410) J,TDI(J),DTD
STOP
XL(J)=T(J)
50 CONTINUE
CALL LINREG(MR)
CALL MEAN
AMU0=EMU(I,1)
CALL PRINT(2)
TMAX=0.
DO 51 J=1,MR
51 TMAX= AMAX1(TMAX,T(J))
LINT=(400-5*MR)/2
LMAX=LINT+3*MR
DO 52 L=1,LINT
X(L)=1.+(TMAX-1.)*(L-1)/FLOAT(LINT-1)
Y(L)=A*X(L)+B
NB(L)=5
52 CONTINUE
DO 53 J=1,MR
K=3*(J-1)+LINT+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
53 CONTINUE
INIT=1
DCHI=(DAF(I,1)/1000.)**2.
DO 55 J=1,MR
JR=J
AMU0=EMU(I,J)
CALL LSAFTAU
55 CONTINUE
INIT=0
AMU0=EMU(I,1)
C
60 CONTINUE
LS=1
ERR=.F.
MU=AMIN1(AMU,AMU0)
DMU=ABS(AMU-AMU0)
61 IF(MU.LE.0.) GOTO 96
62 CALL ROOTS(MU,DMU)
MU=MU-DMU
DMU=3.*DMU
IF(ERRMU) GOTO 61
MU=MU+DMU/3.
MU=FMU(MU)
CALL MEAN
CALL PRINT(3)
DO 63 J=1,MR
K=LMAX+J
X(K)=T(J)
Y(K)=ETA0(J)
NB(K)=3
63 CONTINUE
LMAX=LMAX+MR
C
70 CONTINUE
LS=0
ERR=.F.
MU=AMU
DMUS=SQRT(SA)/14.69
DMU=DMUS/3.
MU=MU-DMU/2.
71 MU=MU-DMU
IF(MU.LE.0.) GOTO 97
72 DMU=3.*DMU
CALL ROOTS(MU,DMU)
IF(ERRMU) GOTO 71
MU=FMU(MU)
CALL PRINT(4)
C
80 CONTINUE
LMIN=LMAX+1
LMAX=LMAX+LINT
DO 81 L=LMIN,LMAX
X(L)=1.+(TMAX-1.)*(L-LMIN)/FLOAT(LINT-1)
Y(L)=A*X(L)+B
NB(L)=6
81 CONTINUE
DO 82 J=1,MR
K=LMAX+J
X(K)=T(J)
Y(K)=ETA0(J)
NB(K)=4
82 CONTINUE
LMAX=LMAX+MR
GOTO 20
C
96 IF(ERR) GOTO 90
MU=EMU(I,1)/10.
ERR=.T.
GOTO 62
97 IF(ERR) GOTO 90
MU=EMU(I,1)/10.
ERR=.T.
GOTO 72
90 CONTINUE
WRITE(4,200)
C
20 CALL GRAPH(X,Y,NB,LMAX)
C
10 CONTINUE
SUBROUTINE INPUT
C
COMMON/INI/ ANF(10,10),DANF(10,10),PH(10,10),DPH(10,10)
$,AF(10,10),DAF(10,10),ETA(10,10),DETA(10,10),EMU(10,10),
$T(10),DT(10),ANFM(10),DANFM(10)
CDC $,F(10,512)
COMMON/TESPAR/ A0M,DA0M,AMU0,AF0(10),ETA0(10),A0(10)
$,TDI(10),TDP(10)
COMMON/PARAM/IDENT,HANG,HI,HF,DK,NP,PI,SQHI,DKHI,I0DH,H0AL
COMMON/STEER/MR,ADENT(10),MF,FREQ(10),I0
COMMON /LINE/ A,B,DA,DB,SA,SB,SAB,TD,AMU,RR,RMSD,CHI
COMMON /REGRES/ XL(10),YL(10),DY(10)
COMMON/FRERUN/ IF,JR,INIT,LS
COMMON/ACCUR/ DDMU1,DDMU2,DCHI,ERRMU
C
REAL APM(10)
LOGICAL ERRMU
610 FORMAT(A2)
620 FORMAT(A6)
630 FORMAT(I2)
640 FORMAT(8(G10.3))
650 FORMAT(A3,A1,A2)
C
PI=4*ATAN(1.)
DDMU1=DDMU2=1.E-4
ICHAN=6
DO 5 J=1,10
5 APM(J)=1.
DO 10 J=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(J)) GO TO 2
BACKSPACE ICHAN
READ(ICHAN,650) IMIST,IDENT,IMIST
READ(ICHAN,630) IMIST
READ(ICHAN,640) DK,HANG,HI,AMIST,T(J),DT(J),AMIST,FNP
IMIST=IMIST-8
IF (IMIST.GT.0) READ(ICHAN,640) (APM(I),I=1,IMIST)
NP=IFIX(FNP)
DK=DK*1.E-3
READ(ICHAN,610) AMIST
READ(ICHAN,630) IMIST
READ(ICHAN,630) NFREQ
DO 20 K=1,NFREQ
READ(ICHAN,640) QN1,QN2,QN3,QN4,DQN1,DQN2,DQN3,DQN4
DO 30 I=1,MF
IF (ABS(QN1-FREQ(I)).GT.1./(DK*NP)) GO TO 30
ANF(I,J)=QN1
PH(I,J)=QN2
AF(I,J)=QN3
ETA(I,J)=QN4
DANF(I,J)=DQN1
DPH(I,J)=DQN2
DAF(I,J)=DQN3
DETA(I,J)=DQN4
EMU(I,J)=APM(K-1)*1.E-3*ANF(I,J)
30 CONTINUE
20 CONTINUE
10 CONTINUE
HF=HI/(1+DK*(NP-1)*HI)
SQHI=SQRT(HI)
DKHI=DK*HI
I0DH=(1.+I0*DKHI)
H0AL=HI/14.69/I0DH
I0DH=SQRT(I0DH)
IF(IDENT.NE.1HM)SQHI=1./SQHI
DO 40 I=1,MF
SA=SB=SP=.0
DO 50 J=1,MR
ANFP=(1./DANF(I,J))**2.
SA=SA+ANF(I,J)*ANFP
SB=SB+ANF(I,J)**2*ANFP
SP=SP+ANFP
50 CONTINUE
ANFM(I)=SA/SP
DANFM(I)=SQRT((SB/SP-ANFM(I)**2)/(MR-1))
40 CONTINUE
RETURN
END
SUBROUTINE MEANTD
C
COMMON/INI/ ANF(10,10),DANF(10,10),PH(10,10),DPH(10,10)
$,AF(10,10),DAF(10,10),ETA(10,10),DETA(10,10),EMU(10,10),
$T(10),DT(10),ANFM(10),DANFM(10)
CDC $,F(10,512)
COMMON/TESPAR/ A0M,DA0M,AMU0,AF0(10),ETA0(10),A0(10)
$,TDI(10),TDP(10)
COMMON/PARAM/IDENT,HANG,HI,HF,DK,NP,PI,SQHI,DKHI,I0DH,H0AL
COMMON/STEER/MR,ADENT(10),MF,FREQ(10),I0
COMMON /LINE/ A,B,DA,DB,SA,SB,SAB,TD,AMU,RR,RMSD,CHI
COMMON /REGRES/ XL(10),YL(10),DY(10)
COMMON/FRERUN/ IF,JR,INIT,LS
COMMON/ACCUR/ DDMU1,DDMU2,DCHI,ERRMU
C
SP=0.
SA=0.
SB=0.
DO 10 J=1,MR
SA=SA+TD(J)*TDP(J)
SB=SB+TD(J)**2*TDP(J)
SP=SP+TDP(J)
10 CONTINUE
TDM=SA/SP
DTDM=SQRT((SB/SP-TDM**2)/(MR-1))
RETURN
END
SUBROUTINE MEAN
C
COMMON/INI/ ANF(10,10),DANF(10,10),PH(10,10),DPH(10,10)
$,AF(10,10),DAF(10,10),ETA(10,10),DETA(10,10),EMU(10,10),
$T(10),DT(10),ANFM(10),DANFM(10)
CDC $,F(10,512)
COMMON/TESPAR/ A0M,DA0M,AMU0,AF0(10),ETA0(10),A0(10)
$,TDI(10),TDP(10)
COMMON/PARAM/IDENT,HANG,HI,HF,DK,NP,PI,SQHI,DKHI,I0DH,H0AL
COMMON/STEER/MR,ADENT(10),MF,FREQ(10),I0
COMMON /LINE/ A,B,DA,DB,SA,SB,SAB,TD,AMU,RR,RMSD,CHI
COMMON /REGRES/ XL(10),YL(10),DY(10)
COMMON/FRERUN/ IF,JR,INIT,LS
COMMON/ACCUR/ DDMU1,DDMU2,DCHI,ERRMU
C
SP=0.
SA=0.
SB=0.
DO 10 J=1,MR
TAU=ETA0(J)*DK*NP
EX=EXP(-TAU)
A0(J)=AF0(J)*TAU*SQHI*EXP(ETA0(J)/HI)/(T(J)*(1.-EX))
PDETA0=(1.-TAU*EX/(1.-EX))/ETA0(J)+1./HI
A0P=A0(J)**2*((DAF(I,J)/AF0(J))**2+(DETA(I,J)*PDETA0)**2)
A0P=1./A0P
SA=SA+A0(J)*A0P
SB=SB+A0(J)**2*A0P
SP=SP+A0P
10 CONTINUE
A0M=SA/SP
DA0M=SQRT((SB/SP-A0M**2)/(MR-1))
RETURN
END
FUNCTION FMU(MU)
C
COMMON/INI/ ANF(10,10),DANF(10,10),PH(10,10),DPH(10,10)
$,AF(10,10),DAF(10,10),ETA(10,10),DETA(10,10),EMU(10,10),
$T(10),DT(10),ANFM(10),DANFM(10)
CDC $,F(10,512)
COMMON/TESPAR/ A0M,DA0M,AMU0,AF0(10),ETA0(10),A0(10)
$,TDI(10),TDP(10)
COMMON/PARAM/IDENT,HANG,HI,HF,DK,NP,PI,SQHI,DKHI,I0DH,H0AL
COMMON/STEER/MR,ADENT(10),MF,FREQ(10),I0
COMMON /LINE/ A,B,DA,DB,SA,SB,SAB,TD,AMU,RR,RMSD,CHI
COMMON /REGRES/ XL(10),YL(10),DY(10)
COMMON/FRERUN/ IF,JR,INIT,LS
COMMON/ACCUR/ DDMU1,DDMU2,DCHI,ERRMU
C
1000 FORMAT(5X'MU ='E17.10' FMU ='E17.10' AMU ='E17.10)
REAL MU
LOGICAL ERRMU
AMU0=MU
DO 10 J=1,MR
JR=J
IF(LS.EQ.0) CALL LSTAU
IF(LS.EQ.1) CALL LSAFTAU
YL(J)=ETA0(J)
DY(J)=DETA(IF,J)
XL(J)=T(J)
10 CONTINUE
CALL LINREG(MR)
FMU=AMU-MU
C WRITE(4,1000) MU,FMU,AMU
RETURN
END
SUBROUTINE LSAFTAU
C
COMMON/INI/ ANF(10,10),DANF(10,10),PH(10,10),DPH(10,10)
$,AF(10,10),DAF(10,10),ETA(10,10),DETA(10,10),EMU(10,10),
$T(10),DT(10),ANFM(10),DANFM(10)
CDC $,F(10,512)
COMMON/TESPAR/ A0M,DA0M,AMU0,AF0(10),ETA0(10),A0(10)
$,TDI(10),TDP(10)
COMMON/PARAM/IDENT,HANG,HI,HF,DK,NP,PI,SQHI,DKHI,I0DH,H0AL
COMMON/STEER/MR,ADENT(10),MF,FREQ(10),I0
COMMON /LINE/ A,B,DA,DB,SA,SB,SAB,TD,AMU,RR,RMSD,CHI
COMMON /REGRES/ XL(10),YL(10),DY(10)
COMMON/FRERUN/ IF,JR,INIT,LS
COMMON/ACCUR/ DDMU1,DDMU2,DCHI,ERRMU
C
1000 FORMAT(5X,A6' CHISQ ='E17.10)
Q3=AF0(JR)
Q4=ETA0(JR)*DK*NP
BE=2.*14.69*AMU0*T(JR)
EBD=EXP(-BE*DK)
H2=1./NP
CHISQ=.0
10 FLD1OV=1.
EBH=EXP(-BE/HI)
E8=E89=1./(1.-EXP(-Q4))
E9=EXP(-Q4/NP)
E3=1.-E8
IF(INIT.EQ.0) GOTO 50
DO 20 K=1,NP
SQTFLD=SQRT(FLD1OV)
IF(IDENT.NE.1HM) SQTFLD=1./SQTFLD
AM=Q3*SQTFLD/(1.-EBH)
EX=Q4*E89
F(JR,K)=AM*EX
FLD1OV=FLD1OV+DKHI
EBH=EBH*EBD
E89=E89*E9
20 CONTINUE
RETURN
50 BR1=BR2=.0
AR11=AR12=AR22=.0
CHISQLA=CHISQ
CHISQ=.0
DO 60 K=1,NP
SQTFLD=SQRT(FLD1OV)
IF(IDENT.NE.1HM) SQTFLD=1./SQTFLD
AM=Q3*SQTFLD/(1.-EBH)
EX=Q4*E89
R1=F(JR,K)-AM*EX
CHISQ=CHISQ+R1**2
DFDQ3=AM*EX/Q3
DFDQ4=AM*E89*(1.+Q4*E3)
BR1=BR1+R1*DFDQ3
BR2=BR2+R1*DFDQ4
AR11=AR11+DFDQ3**2
AR22=AR22+DFDQ4**2
AR12=AR12+DFDQ3*DFDQ4
FLD1OV=FLD1OV+DKHI
EBH=EBH*EBD
E89=E89*E9
E3=E3-H2
60 CONTINUE
CHISQ=CHISQ/(NP-2)
C WRITE(4,1000) ADENT(JR),CHISQ
IF(ABS(CHISQ-CHISQLA).LT.DCHI) GOTO 100
DET=AR11*AR22-AR12**2
Q3=Q3+(BR1*AR22-BR2*AR12)/DET
Q4=Q4+(BR2*AR11-BR1*AR12)/DET
GOTO 10
100 AF0(JR)=Q3
ETA0(JR)=Q4/(DK*NP)
RETURN
END
SUBROUTINE LSTAU
C
COMMON/INI/ ANF(10,10),DANF(10,10),PH(10,10),DPH(10,10)
$,AF(10,10),DAF(10,10),ETA(10,10),DETA(10,10),EMU(10,10),
$T(10),DT(10),ANFM(10),DANFM(10)
CDC $,F(10,512)
COMMON/TESPAR/ A0M,DA0M,AMU0,AF0(10),ETA0(10),A0(10)
$,TDI(10),TDP(10)
COMMON/PARAM/IDENT,HANG,HI,HF,DK,NP,PI,SQHI,DKHI,I0DH,H0AL
COMMON/STEER/MR,ADENT(10),MF,FREQ(10),I0
COMMON /LINE/ A,B,DA,DB,SA,SB,SAB,TD,AMU,RR,RMSD,CHI
COMMON /REGRES/ XL(10),YL(10),DY(10)
COMMON/FRERUN/ IF,JR,INIT,LS
COMMON/ACCUR/ DDMU1,DDMU2,DCHI,ERRMU
C
1000 FORMAT(5X,A6' CHISQ ='E17.10)
Q4=ETA0(JR)*DK*NP
BE=2.*14.69*AMU0*T(JR)
EBD=EXP(-BE*DK)
H2=1./NP
CHISQ=.0
10 FLD1OV=1.
EBH=EXP(-BE/HI)
E9=EXP(-Q4/NP)
EX=EXP(-Q4/(HI*DK*NP))
H3=1./(HI*DK*NP)
BR=.0
AR=.0
CHISQLA=CHISQ
CHISQ=.0
DO 20 K=1,NP
SQTFLD=SQRT(FLD1OV)
IF(IDENT.NE.1HM) SQTFLD=1./SQTFLD
AM=A0M*T(JR)*SQTFLD/(SQHI*(1.-EBH))
R1=F(JR,K)-AM*EX
CHISQ=CHISQ+R1**2
DFDQ4=-H3*AM*EX
BR=BR+R1*DFDQ4
AR=AR+DFDQ4*DFDQ4
FLD1OV=FLD1OV+DKHI
EBH=EBH*EBD
EX=EX*E9
H3=H3+H2
20 CONTINUE
CHISQ=CHISQ/(NP-1)
C WRITE(4,1000) ADENT(JR),CHISQ
IF(ABS(CHISQ-CHISQLA).LT.DCHI) GOTO 30
Q4=Q4+BR/AR
GOTO 10
30 ETA0(JR)=Q4/(DK*NP)
RETURN
END
SUBROUTINE PRINT(IFUN)
C
COMMON/INI/ ANF(10,10),DANF(10,10),PH(10,10),DPH(10,10)
$,AF(10,10),DAF(10,10),ETA(10,10),DETA(10,10),EMU(10,10),
$T(10),DT(10),ANFM(10),DANFM(10)
CDC $,F(10,512)
COMMON/TESPAR/ A0M,DA0M,AMU0,AF0(10),ETA0(10),A0(10)
$,TDI(10),TDP(10)
COMMON/PARAM/IDENT,HANG,HI,HF,DK,NP,PI,SQHI,DKHI,I0DH,H0AL
COMMON/STEER/MR,ADENT(10),MF,FREQ(10),I0
COMMON /LINE/ A,B,DA,DB,SA,SB,SAB,TD,AMU,RR,RMSD,CHI
COMMON /REGRES/ XL(10),YL(10),DY(10)
COMMON/FRERUN/ IF,JR,INIT,LS
COMMON/ACCUR/ DDMU1,DDMU2,DCHI,ERRMU
C
LOGICAL ERRMU
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 : ',F10.3,' T',T50'MU*P='
$,T60,3PF7.3,T70,'1E-3*1/T')
540 FORMAT(//2X'RUN'T16'T[K]'T32'FREQ[T]'T51'PHASE'T71'AMP'T91,
$'ETA[T]'T107'A0'/)
550 FORMAT(1X,A6,2X,0PF6.3'+/-'F5.3,2X,F9.2'+/-'F5.2,3X,F6.1,
$'+/-'F4.1,3X,1PE9.2'+/-'E9.2,3X,0PF6.2'+/-'F5.2,3X,1PE9.2)
560 FORMAT(1X,//,T8,'MEAN FREQUENCY[T]='F9.2'+/-'F5.2,T95,
$'MEAN A0='1PE9.2,/,T8,'PERIOD[1/T]= '1PE9.2,
$T95'VAR.A0= '1PE9.2,////)
570 FORMAT(1X,T8,'ETA FIT:',T20,'X =',T30,F7.3,T40,'K',T70,
$'MU =',T80,F7.3,/,T8,'........',/,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,
$'CHI=',T80,0PF7.3,/,T70,'RMSD=',T80,F7.3,/,T70,
$'RR=',T80,F7.3,//)
590 FORMAT(1X,T8,'ETA+ FIT:',T20,'X =',T30,F7.3,T40,'K',T70,
$'MU =',T80,F7.3,//,T33,'+/-',T38,F7.3,
$T83,'+/-',T88,F7.3,' (LINE FIT)',/,T38,F7.3,T88,F7.3,
$' (VARIANCE)',,//)
600 FORMAT(1X,T8,'ETA0 FIT:',T20,'X =',T30,F7.3,T40,'K',T70,
$'MU =',T80,F7.3,/,T8,'*********',/,T33,'+/-',T38,F7.3,
$T83,'+/-',T88,F7.3,' (LINE FIT)',/,T38,F7.3,T88,F7.3,
$' (VARIANCE)',,//)
C
C !!!!!!!!!!!!!!!!!!!!!!!!!
IVERS=8
C !!!!!!!!!!!!!!!!!!!!!!!!!
IF (IFUN.EQ.2) GO TO 2
IF (IFUN.EQ.3) GO TO 3
IF (IFUN.EQ.4) GO TO 4
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)
3 AMP=AMU0/ANFM(I)
WRITE(4,530) FREQ(I),AMP
GOTO 6
4 DA0M=.0
DO 5 J=1,MR
A0(J)=A0M
TAU=ETA0(J)*DK*NP
EX=EXP(-TAU)
AF0(J)=A0M*T(J)*(1.-EX)/(TAU*SQHI*EXP(ETA0(J)/HI))
5 CONTINUE
6 WRITE(4,540)
WRITE(4,550) (ADENT(J),T(J),DT(J),ANF(I,J),DANF(I,J),
1PH(I,J),DPH(I,J),AF0(J),DAF(I,J),ETA0(J),DETA(I,J),A0(J),
2J=1,MR)
PER=1./ANFM(I)
WRITE(4,560) ANFM(I),DANFM(I),A0M,PER,DA0M
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,570) TD,AMU,DTDL,DMUL,DTDS,DMUS
IF (IFUN.EQ.3) WRITE(4,590) TD,AMU,DTDL,DMUL,DTDS,DMUS
IF (IFUN.EQ.4) WRITE(4,600) TD,AMU,DTDL,DMUL,DTDS,DMUS
AMP=AMU/ANFM(I)
WRITE(4,580) AMP,CHI,RMSD,RR
IF (IFUN.GE.3) WRITE(4,500)
RETURN
END
SUBROUTINE ROOTS(ROOT,DR)
C
COMMON/INI/ ANF(10,10),DANF(10,10),PH(10,10),DPH(10,10)
$,AF(10,10),DAF(10,10),ETA(10,10),DETA(10,10),EMU(10,10),
$T(10),DT(10),ANFM(10),DANFM(10)
CDC $,F(10,512)
COMMON/TESPAR/ A0M,DA0M,AMU0,AF0(10),ETA0(10),A0(10)
$,TDI(10),TDP(10)
COMMON/PARAM/IDENT,HANG,HI,HF,DK,NP,PI,SQHI,DKHI,I0DH,H0AL
COMMON/STEER/MR,ADENT(10),MF,FREQ(10),I0
COMMON /LINE/ A,B,DA,DB,SA,SB,SAB,TD,AMU,RR,RMSD,CHI
COMMON /REGRES/ XL(10),YL(10),DY(10)
COMMON/FRERUN/ IF,JR,INIT,LS
COMMON/ACCUR/ DDMU1,DDMU2,DCHI,ERRMU
C
LOGICAL FEHLER
FEHLER=.F.
ZA=ROOT
ZB=ROOT+DR
FA=FMU(ZA)
SIG=SIGN(1.,FA)
FA=SIG*FA
FB=SIG*FMU(ZB)
IF(FA*FB) 30,30,500
C
C LINEARE INTERPOLATION UND INTERVALLSCHACHTELUNG
C
30 ZC=(ZA*FB-ZB*FA)/(FB-FA)
FC=SIG*FMU(ZC)
IF(ABS(FC).LT.DDET) GOTO 91
IF(FC.LT.0.) GOTO 70
ZD=(ZB+ZC)/2.
FD=SIG*FMU(ZD)
50 IF(ABS(FD).LT.DDET) GOTO 92
IF(FD.LT.0.) GOTO 100
ZC=ZD $ FC=FD
ZD=(ZB+ZD)/2.
FD=SIG*FMU(ZD)
GOTO 50
70 ZD=ZC $ FD=FC
ZC=(ZA+ZC)/2.
FC=SIG*FMU(ZC)
IF(ABS(FC).LT.DDET) GOTO 91
IF(FC.LT.0.) GOTO 70
100 ZA=ZC $ ZB=ZD
FA=FC $ FB=FD
GOTO 96
91 ZA=ZB=ZC
FA=FB=.0
GOTO 96
92 ZA=ZB=ZD
FA=FB=.0
96 CONTINUE
IF(ABS(FA-FB).GT.DDET.AND.ABS(ZA-ZB).GT.DK) GOTO 30
ROOT=(ZA+ZB)/2.
RETURN
500 FEHLER=.T.
RETURN
END
SUBROUTINE LINREG(MR)
C
COMMON/INI/ ANF(10,10),DANF(10,10),PH(10,10),DPH(10,10)
$,AF(10,10),DAF(10,10),ETA(10,10),DETA(10,10),EMU(10,10),
$T(10),DT(10),ANFM(10),DANFM(10)
CDC $,F(10,512)
COMMON/TESPAR/ A0M,DA0M,AMU0,AF0(10),ETA0(10),A0(10)
$,TDI(10),TDP(10)
COMMON/PARAM/IDENT,HANG,HI,HF,DK,NP,PI,SQHI,DKHI,I0DH,H0AL
COMMON/STEER/MR,ADENT(10),MF,FREQ(10),I0
COMMON /LINE/ A,B,DA,DB,SA,SB,SAB,TD,AMU,RR,RMSD,CHI
COMMON /REGRES/ XL(10),YL(10),DY(10)
COMMON/FRERUN/ IF,JR,INIT,LS
COMMON/ACCUR/ DDMU1,DDMU2,DCHI,ERRMU
C
DIMENSION V(10),P(10)
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
N=MR
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+XL(J)*P(J)
SY=SY+YL(J)*P(J)
SXX=SXX+XL(J)*XL(J)*P(J)
SXY=SXY+XL(J)*YL(J)*P(J)
SYY=SYY+YL(J)*YL(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
TD=B/A
AMU=A/14.69
C VARIANCES FROM LINE FIT
IF (N.GT.2) GO TO 75
DA=0.
DB=0.
RR=0.
RMSD=0.
CHI=0.
GO TO 76
75 CONTINUE
Q=0.
QA=0.
DO 62 J=1,N
V(J)=YL(J)-XL(J)*A-B
Q=Q+V(J)*V(J)*P(J)
QA=QA+(V(J))**2
62 CONTINUE
CHI=SQRT(Q/(N-2))
C=SP*SYY-SY*SY
DA=CHI*SQRT(SP/D)
DB=CHI*SQRT(SXX/D)
RR=SQRT(A*A*D/C)
RMSD=SQRT(QA/(N-2))
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)=1H+
GOTO 60
53 IF(NB(N).NE.4) GOTO 54
LINEA(K)=1H0
GOTO 60
54 IF(NB(N).NE.5) GOTO 55
IF(LINEA(K).NE.1H ) GOTO 60
LINEA(K)=1H.
GOTO 60
55 IF(NB(N).NE.6) GOTO 56
IF(LINEA(K).NE.1H ) GOTO 60
LINEA(K)=1H*
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