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