C ************************************** C Main Programm peak.f * C ************************************** C Yui IGARASHI & Kenji SUZUKI CHARACTER*18 FNAME REAL*8 A,B,X,Y,AY,BY,H,RATIO,XA,XB,XP,YP REAL*8 PEAK,PRT,TEMPX,TEMPY,Z INTEGER*2 N,I DIMENSION X(1:1000),Y(1:1000),TEMPX(1:1000),TEMPY(1:1000) DIMENSION Z(1:11) WRITE(6,*) 'INPUT DATA FILE NAME ' READ(5,*) FNAME OPEN(UNIT=2,FILE=FNAME,STATUS='OLD') C _______ PRT: PRESET TIME (SEC) READ(2,*) PRT I=1 C ____Y:X-RAY COUNTS, X:2-THETA (DEG) FROM HIGH ANGLE TO LOW ANGLE 200 READ(2,*,END=201) Z(1),Z(2),Z(3),Z(4),Z(5),Z(6),Z(7), 1 Z(8),Z(9),Z(10),Z(11) IF (I .EQ. 1) THEN ORG=Z(10) END IF X(I)=Z(1) Y(I)=Z(11)*ORG/Z(10)/Z(7) WRITE(6,*) X(I),Y(I) I=I+1 GOTO 200 201 CLOSE(UNIT=2) N=I-1 C ____X(I)=2TH, Y(I) (cps) 667 CONTINUE CALL BACKG(X,Y,N,A,B) C WRITE(6,*) 'A= ',A,' B= ',B C __________PREPROCESSED CURVE______ C ____GNUPLOT 1______________ OPEN(UNIT=2,FILE='bg.gnu',STATUS='UNKNOWN') DO 202 I=1,N WRITE(2,*) X(I),Y(I) 202 CONTINUE WRITE(2,*) WRITE(2,*) AY=A*X(1)+B BY=A*X(N)+B WRITE(2,*) X(1),AY WRITE(2,*) X(N),BY CLOSE(UNIT=2) OPEN(UNIT=2,FILE='x-data.frn',STATUS='UNKNOWN') C DO 10 I=1,N Y(I)=Y(I)-A*X(I)-B WRITE(2,100) X(I),Y(I) 10 CONTINUE CLOSE(UNIT=2) 100 FORMAT(F10.4,2X,F10.3) C ________ PEAK TOP FITTED BY PARABOLA __________ C PEAK POSITION (XP,YP) CALL TOP(X,Y,N,XP,YP) WRITE(6,*) 'INPUT RATIO TO PEAK HEIGHT ' READ(5,*) RATIO H=RATIO*YP C _____ BOTH SLOPES APPROXIMATED BY CUBIC CALL CUBIC(X,Y,N,H,XA,XB) WRITE(6,*) XA,XB PEAK=(XA+XB)/2.D0 WRITE(6,*) '2_theta (deg)=',PEAK WRITE(6,*) 'Width (deg)=',XB-XA C ____GNUPLOT 2______________ OPEN(UNIT=2,FILE='peak.gnu',STATUS='UNKNOWN') DO 203 I=1,N WRITE(2,*) X(I),Y(I) 203 CONTINUE WRITE(2,*) WRITE(2,*) WRITE(2,*) XA,H WRITE(2,*) XB,H WRITE(2,*) WRITE(2,*) WRITE(2,*) PEAK,YP WRITE(2,*) PEAK,0.0 CLOSE(UNIT=2) STOP END C _/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/ C BACK GROUND C ______________________________ SUBROUTINE BACKG(X,Y,N,A,B) C REAL*8 A,B,X,Y,C,MA,MB INTEGER*2 N,I,NM,NN DIMENSION X(1:1000),Y(1:1000) DIMENSION MA(1:100,1:100),MB(1:100) NM=5 NN=2 C C=0.D0 DO 30 I=1,NM C=C+X(I)**2 30 CONTINUE DO 31 I=N-NM+1,N C=C+X(I)**2 31 CONTINUE MA(1,1)=C C=0.D0 DO 40 I=1,NM C=C+X(I) 40 CONTINUE DO 41 I=N-NM+1,N C=C+X(I) 41 CONTINUE MA(1,2)=C MA(2,1)=MA(1,2) MA(2,2)=FLOAT(NM*2) C C=0.D0 DO 20 I=1,NM C=C+X(I)*Y(I) 20 CONTINUE DO 21 I=N-NM+1,N C=C+X(I)*Y(I) 21 CONTINUE MB(1)=C C=0.D0 DO 50 I=1,NM C=C+Y(I) 50 CONTINUE DO 51 I=N-NM+1,N C=C+Y(I) 51 CONTINUE MB(2)=C C CALL GAUSS(MA,MB,NN) A=MB(1) B=MB(2) RETURN END C___________________________________________________ * Solving simultaneous equations * using Gauss method * CALL GAUSS C___________________________________________________ SUBROUTINE GAUSS(A,B,N) REAL*8 A,B,P,Q INTEGER*2 I,J,K,L,N DIMENSION A(1:100,1:100),B(1:100) C FRONT PROCESS C SETTING L-LINE DO 10 L=1,N P=A(L,L) DO 20 I=L,N A(L,I)=A(L,I)/P 20 CONTINUE B(L)=B(L)/P C SUBSTITUTION DO 30 I=L+1,N Q=A(I,L) DO 40 J=L,N A(I,J)=A(I,J)-Q*A(L,J) 40 CONTINUE B(I)=B(I)-Q*B(L) 30 CONTINUE 10 CONTINUE C POST PROCESS C SETTING L-LINE DO 100 L=N,1,-1 P=A(L,L) DO 110 I=N,L,-1 A(L,I)=A(L,I)/P 110 CONTINUE B(L)=B(L)/P C SUBSTITUTION DO 120 I=L-1,1,-1 Q=A(I,L) DO 130 J=N,I,-1 A(I,J)=A(I,J)-Q*A(L,J) 130 CONTINUE B(I)=B(I)-Q*B(L) 120 CONTINUE 100 CONTINUE RETURN END C _____________________________________________ C C Peak top approximated by parabola SUBROUTINE TOP(X,Y,N,X0,Y0) REAL*8 X,Y,Ymax,M,D,X0,Y0,SUM REAL*8 X1,Y1,A,B,C REAL*8 C1,C2,C3,C4,C5,C6,C7 INTEGER*2 N,I,V,J,Imax DIMENSION X(1:1000),Y(1:1000),M(1:100,1:100),D(1:100) C _____ V: NUMBER OF UNKNOWN VARIABLES ___________ V=3 J=7 OPEN(UNIT=1,FILE='x-data.frn',STATUS='UNKNOWN') C I=1 100 READ(1,*,END=101) X(I),Y(I) I=I+1 GOTO 100 101 CLOSE(UNIT=1) N=I-1 C Ymax=Y(1) Imax=1 DO 200 I=2,N IF (Y(I).GE.Ymax) THEN Ymax=Y(I) Imax=I END IF 200 CONTINUE C WRITE(6,*) Imax,X(Imax),Ymax C C1=0.D0 C2=0.D0 C3=0.D0 C4=0.D0 C5=0.D0 C6=0.D0 C7=0.D0 DO 10 I=Imax-3,Imax+3 C1=C1+X(I)**4 C2=C2+X(I)**3 C3=C3+X(I)**2 C4=C4+X(I) C5=C5+Y(I)*X(I)**2 C6=C6+Y(I)*X(I) C7=C7+Y(I) 10 CONTINUE M(1,1)=C1 M(1,2)=C2 M(1,3)=C3 M(2,1)=M(1,2) M(2,2)=M(1,3) M(2,3)=C4 M(3,1)=M(1,3) M(3,2)=M(2,3) M(3,3)=FLOAT(J) C D(1)=C5 D(2)=C6 D(3)=C7 C CALL GAUSS(M,D,V) C A=D(1) B=D(2) C=D(3) X0=-B/2.D0/A Y0=(4.D0*A*C-B**2)/4.D0/A WRITE(6,*) 'PEAK TOP= ',X0,Y0 C ____________GNUPLOT OPEN(UNIT=2,FILE='top.gnu',STATUS='UNKNOWN') C DO 900 I=1,N WRITE(2,*) X(I),Y(I) 900 CONTINUE WRITE(2,*) WRITE(2,*) DO 990 X1=X0-0.5,X0+0.5,0.01 Y1=A*X1**2+B*X1+C IF (Y1 .GT. 0.E0) THEN WRITE(2,*) X1,Y1 END IF 990 CONTINUE CLOSE(UNIT=2) RETURN END C____________________________________________________ C SLOPE APPROXIMATED BY CUBIC FUNCTION SUBROUTINE CUBIC(X,Y,N,H,XA,XB) REAL*8 X,Y,H,XA,XB,A,B,C,D,MA,MB,XI REAL*8 C1,C2,C3,C4,C5,C6,C7,C8,C9,C10 INTEGER*2 N,I,IA,IB,V DIMENSION X(1:1000),Y(1:1000),XI(0:5) DIMENSION MA(1:100,1:100),MB(1:100) C================== LEFT SIDE SLOPE ============= C IA=1 DO 10 I=2,N IF (Y(I) .GT. H) GOTO 12 10 CONTINUE 12 IA=I write(6,*) IA,X(IA) C1=0.D0 C2=0.D0 C3=0.D0 C4=0.D0 C5=0.D0 C6=0.D0 C7=0.D0 C8=0.D0 C9=0.D0 C10=0.D0 DO 14 I=IA-3,IA+2 c DO 14 I=0,5 C1=C1+X(I)**6 C2=C2+X(I)**5 C3=C3+X(I)**4 C4=C4+X(I)**3 C5=C5+X(I)**2 C6=C6+X(I) C7=C7+Y(I)*X(I)**3 C8=C8+Y(I)*X(I)**2 C9=C9+Y(I)*X(I) C10=C10+Y(I) 14 CONTINUE MA(1,1)=C1 MA(1,2)=C2 MA(1,3)=C3 MA(1,4)=C4 MA(2,1)=MA(1,2) MA(2,2)=MA(1,3) MA(2,3)=MA(1,4) MA(2,4)=C5 MA(3,1)=MA(1,3) MA(3,2)=MA(1,4) MA(3,3)=MA(2,4) MA(3,4)=C6 MA(4,1)=MA(1,4) MA(4,2)=MA(2,4) MA(4,3)=MA(3,4) MA(4,4)=FLOAT(6) MB(1)=C7 MB(2)=C8 MB(3)=C9 MB(4)=C10 V=4 C____________________ GAUSS METHOD CALL GAUSS(MA,MB,V) A=MB(1) B=MB(2) C=MB(3) D=MB(4) write(6,*) A,B,C,D C_______ INTERCEPT POINT CALCULATED BY NEWTON METHOD C XA: LEFT SIDE POINT D=D-H XA=X(IA) CALL NEWTON(A,B,C,D,XA) C================== RIGHT SIDE SLOPE ============= IB=N DO 20 I=N,1,-1 IF (Y(I) .GT. H) GOTO 22 20 CONTINUE 22 IB=I write(6,*) IB,X(IB) C1=0.D0 C2=0.D0 C3=0.D0 C4=0.D0 C5=0.D0 C6=0.D0 C7=0.D0 C8=0.D0 C9=0.D0 C10=0.D0 DO 16 I=IB-2,IB+3 C1=C1+X(I)**6 C2=C2+X(I)**5 C3=C3+X(I)**4 C4=C4+X(I)**3 C5=C5+X(I)**2 C6=C6+X(I) C7=C7+Y(I)*X(I)**3 C8=C8+Y(I)*X(I)**2 C9=C9+Y(I)*X(I) C10=C10+Y(I) 16 CONTINUE MA(1,1)=C1 MA(1,2)=C2 MA(1,3)=C3 MA(1,4)=C4 MA(2,1)=MA(1,2) MA(2,2)=MA(1,3) MA(2,3)=MA(1,4) MA(2,4)=C5 MA(3,1)=MA(1,3) MA(3,2)=MA(1,4) MA(3,3)=MA(2,4) MA(3,4)=C6 MA(4,1)=MA(1,4) MA(4,2)=MA(2,4) MA(4,3)=MA(3,4) MA(4,4)=FLOAT(6) MB(1)=C7 MB(2)=C8 MB(3)=C9 MB(4)=C10 V=4 C____________________ GAUSS METHOD CALL GAUSS(MA,MB,V) A=MB(1) B=MB(2) C=MB(3) D=MB(4) write(6,*) A,B,C,D C_______ INTERCEPT POINT CALCULATED BY NEWTON METHOD C XA: LEFT SIDE POINT D=D-H XB=X(IB) CALL NEWTON(A,B,C,D,XB) RETURN END C__________________NEWTON METHOD SUBROUTINE NEWTON(A,B,C,D,X0) REAL*8 A,B,C,D,X0,X,GAMMA,BETA integer*4 I I=0 10 X=X0 GAMMA=3.D0*A*X*X+2.D0*B*X+C BETA=-2.D0*A*X**3-B*X*X+D X=-BETA/GAMMA IF(ABS(X-X0) .LT. 1.D-4) GOTO 20 X0=X I=I+1 IF(I .GT. 1000) GOTO 20 GO TO 10 20 X0=X RETURN END