Fita uma parábola
program PBFT
C PARABOLIC FITTING PROGRAM
C Y=A+BX+CX**2
C WEIGHTING:W(I)=1/SIGY(I)**2 ENTER SIGY(I)
C X IS THE INDEPENDENT VARIABLE ARRAY WITH N POINTS:INPUT DATA
C Y IS THE DEPENDENT VARIABLE ARRAY WITH N POINTS:INPUT DATA
c w is the weight array
C N IS THE NUMBER OF POINTS: INPUT DATA
C R IS A three LINES ARRAY WHERE:R(1)=A,R(2)=B, R(3)=C:OUTPUT DATA
C SIG2 IS THE VARIANCE OF THE FIT:SIG2=SUM((Y OBS-Y CALC)**2)*P/(N-2):
C OUTPUT DATA
C CM IS THE 3x3 COVARIANCE MATRIX:OUTPUT DATA
C************************************************************************
C Written in 6 jun 86 S.O. Kepler
C Last modified ** 6 jun 86 **
C************************************************************************
implicit double precision(a-h,p-z)
dimension x(1000),y(1000),w(1000),sigy(1000),ycalc(1000)
dimension R(3),A(3,3),B(3),CM(3,3),sigma(3)
DIMENSION P(20),Q(20),RR(20)
CHARACTER*20 infile
321 write(*,322)
322 format(1x,'input file is: ', $)
read(*,323) infile
323 format(a20)
open(2,file=infile,status='old')
c ***********************************
c read in the number of data points
read(2,*)N
c ***********************************
iii=1
write(*,521) iii
521 format(1x,'program completed phase: ', I3)
write(*,522) N
522 format(1x,'N= ',i5)
10 DO 28 I=1,3
DO 78 J=1,3
A(I,J)=0.0d0
B(I)=0.0d0
R(I)=0.0d0
CM(I,J)=0.0d0
78 CONTINUE
28 CONTINUE
c read header***********************************
c *****************************************
c read in data - time (x), intensity (y) and sigma (sigy)
DO 325 I=1,N
325 READ(2,*) X(I),Y(I),sigy(i)
c *****************************************
iii=2
write(*,521) iii
DO 77 I=1,N
XI=X(I)
YI=Y(I)
W(I)=1/SIGY(I)**2
WI=W(I)
A(1,2)=A(1,2)+WI*XI
A(1,3)=A(1,3)+WI*XI**2
A(2,3)=A(2,3)+WI*XI**3
A(3,3)=A(3,3)+WI*XI**4
B(1)=B(1)+WI*YI
B(2)=B(2)+WI*YI*XI
B(3)=B(3)+WI*YI*XI**2
A(1,1)=A(1,1)+WI
77 CONTINUE
A(2,1)=A(1,2)
A(3,1)=A(1,3)
A(2,2)=A(1,3)
A(3,2)=A(2,3)
NN=3
C CALL SYMIN(A,NN)
C SUBROUTINE SYMIN (A,NN)
C
C THIS SUBROUTINE INVERTS A SYMMETRIC MATRIX.
C ***** INPUT DATA *****
C A = A SYMMETRIC MATRIX
C NN = THE RANK OF THE MATRIX. N MUST BE LESS THAN 20.
C ***** OUTPUT DATA *****
C A = THE INVERSE OF THE INPUT MATRIX
C
C DIMENSION A(NN,NN)
C DIMENSION P(20),Q(20),R(20)
1300 FORMAT(13H0SYMIN FAILED)
ZERO = 0.0d0
ONE = 1.0d0
DO 1000 M=1,NN
1000 RR(M) = ONE
DO 200 M=1,NN
BIG = ZERO
DO 3000 L=1,NN
AB = DABS(A(L,L))
IF(AB-BIG)3000,3000,4000
4000 IF(RR(L))1400,3000,1400
1400 BIG = AB
K = L
3000 CONTINUE
IF(BIG)6000,5000,6000
5000 WRITE (*,1300)
GO TO 101
C RETURN
6000 RR(K) = ZERO
Q(K) = ONE/A(K,K)
P(K) = ONE
A(K,K) = ZERO
KM1 = K-1
IF(KM1.EQ.0) GO TO 1600
DO 7000 L=1,KM1
P(L) = A(L,K)
IF(RR(L))9000,8000,9000
8000 Q(L) = A(L,K)*Q(K)
GO TO 7000
9000 Q(L) = -A(L,K)*Q(K)
7000 A(L,K) = ZERO
1600 CONTINUE
KP1 = K+1
IF(KP1.GT.NN) GO TO 1700
DO 1500 L=KP1,NN
IF(RR(L))1200,1100,1200
1200 P(L) = A(K,L)
GO TO 10000
1100 P(L) = -A(K,L)
10000 Q(L) = (-A(K,L))*Q(K)
1500 A(K,L) = ZERO
1700 CONTINUE
DO 200 L=1,NN
DO 200 K=L,NN
200 A(L,K) = A(L,K) + P(L)*Q(K)
M = NN+1
L = NN
DO 27 K=2,NN
M = M-1
L = L-1
DO 27 J=1,L
27 A(M,J) = A(J,M)
C RETURN
C NOW THE A MATRIX IS THE INVERSE OF THE INITIAL ONE
C FINDING R(I),WHERE Y=R(1)+R(2)X+R(3)X**2
DO 199 I=1,3
DO 299 J=1,3
R(I)=R(I)+A(I,J)*B(J)
299 CONTINUE
199 CONTINUE
DO 1 I=1,N
YCALC(I)=R(1)+R(2)*X(I)+R(3)*X(I)**2
1 CONTINUE
C FINDING SIG2
38 CONTINUE
RSIG=0.0d0
DO 1333 I=1,N
RSIG=RSIG+((Y(I)-R(1)-R(2)*X(I)-R(3)*X(I)**2)**2)*W(I)
1333 CONTINUE
SIG2=RSIG/(N-3)
SIG2NEW=RSIG/PTOTAL*N/(N-3)
DO 1001 I=1,3
DO 1002 J=1,3
CM(I,J)=A(I,J)*SIG2NEW
1002 CONTINUE
1001 CONTINUE
do 3915 i=1,3
sigma(i)=dsqrt(cm(i,i))
3915 continue
100 format('Fitting a line y=a+bx+cx**2'/,'a = ',1pe12.5,/'b = ',1pe
&12.5,/'c = ',1pe12.5,/'sig2= ',1pe12.3)
c... output
c
c write header*******************************
c
c *******************************************
432 write(*,433)
433 format(1x,'output file name is: ', $)
read(*,323) infile
open(3,file=infile,status='new')
write(3,99)N
99 format(' npoints= ',i5)
c output results
write(3,100)r,sig2
write(3,110)(sigma(i),i=1,3)
110 format('siga =',1pe12.3,/'sigb =',1pe12.3,/'sigc =',1pe12.3)
write(3,107)((cm(i,j),j=1,3),i=1,3)
107 format(31x,'Correlation Matrix'/,
&29x,'a, b, c'/,
&1x,3(3(0pg18.6,2x)/,1x))
101 continue
c **********************
stop
end