c
c Programa que calcula un ajuste gaussiano para el espectro de
c Para ello se utiliza el archivo resultado de la sustración de fondo del c programa
c difesp.f, del que obtenemos un archivo compuesto por:
c -Primera fila: tiempo de medida del espectro
c -Segunda fila: numero de datos del archiv
c -Por último, número de cuentas en cada canal y desviación
c
c character*9 me indica que los archivos que vienen a continuación
c deban tener 9 caracteres (sin contar .ext), mientras que ext1*4 implica c que debe
c tener una extensión de longitud 4 incluido el .
c A continuación dimensionamos todas las variables, para ello le damos c dimensión
c a nmax, y todo lo demas va con su dimensión
c
character*9 name1,ext1*4,ext2*4
c
parameter(nmax=1200)
dimension ydat(nmax),pac(3),lista(3),covar(3,3),afun(10)
c
c
c Pedimos por pantalla los espectros que deseamos manipular
c
write(*,*)'Programa que realiza un ajuste gaussiano'
write(*,*)'de un espectro quantulus'
write(*,*)'Nombre fichero 1'
read(*,*)name1
c
c
c Definimos que extensión corresponde a ext1 y a ext2
c
ext1='.dat'
ext2='.res'
c
c
c
c Abrimos el archivo
c
open(unit=1,file=name1//ext1)
open(unit=2,file=name1//ext2)
c
c Le pedimos que lea del archivo de entrada el tiempo de medida y el número
c dedatos, para que a continuación lea los datos por canal con un bucle que
c lea desde el canal 1 hasta el impuesto en el archivo como ndata
c
read(1,*)t
write(*,*)t
read(1,*)ndat
write(*,*)ndat
do 11 n=1,ndat
read(1,*)ydat(n),sydat(n)
write(*,*)n,ydat(n),sydat(n)
11 continue
c
c
c
c Le pedimos que realice el ajuste gaussiano
c Datos relacionados con el ajuste
c
ma=3
mfit=ma
ncvm=3
ma=1
mfit=ma
lista(1)=1
lista(2)=2
lista(3)=3
c
call lfit(xf,ydat,sigydat,ndata,af,ma,lista,covar,mfit,ncvm,chisq)
c
c
write(2,*)'AJUSTE GAUSIANO'
write(2,*)'ESPECTRO QUANTULUS'
write(2,*)(afun(m),m=1,ma)
write(2,*)chisq/(lm-1)
c
c
c Calculo de A,B y C
c
A=afun(1)
sigA=sqrt(covar(1,1))*A/afun(1)
c
write(2,*)A,sigA
c
B=afun(2)
sigA=sqrt(covar(2,2))*B/afun(2)
c
write(2,*)B,sigB
c
c=sqrt(afun(3)/2)
sigA=sqrt(covar(3,3))*B/afun(3)
c
write(2,*)C,sigC
c
write(18,*)chisq/(ndata-1)
write(2,*)'Resultados: chisqr,A,sA,B,sB,C,sC,covar(1,2)'
do 14 n=1,ndata
call funcs(ydat(n),afun,ma)
ydat(n)=afun(1)*exp(((n-afun(2))/2/afun(3))**2)
write(18,*)chisqr,A,sA,B,sB,C,sC,covar(1,2)
14 continue
c
c
end
c
SUBROUTINE LFIT(xf,YF,SIG,NF,Af,MA,LISTA,COVAR,MFIT,NCVM,CHISQ)
c
DIMENSION xF(NF),YF(NF),SIG(NF),Af(MA),LISTA(MFIT),
* COVAR(NCVM,NCVM),BETA(100),AFUNC(100)
KK=MFIT+1
DO 12 J=1,MA
IHIT=0
DO 11 K=1,MFIT
IF (LISTA(K).EQ.J) IHIT=IHIT+1
11 CONTINUE
IF (IHIT.EQ.0) THEN
LISTA(KK)=J
KK=KK+1
ELSE IF (IHIT.GT.1) THEN
PAUSE 'Improper set in LISTA'
ENDIF
12 CONTINUE
IF (KK.NE.(MA+1)) PAUSE 'Improper set in LISTA'
DO 14 J=1,MFIT
DO 13 K=1,MFIT
COVAR(J,K)=0.
13 CONTINUE
BETA(J)=0.
14 CONTINUE
DO 18 I=1,nf
CALL FUNCS(xF(I),AFUNC,MA)
YM=YF(I)
IF(MFIT.LT.MA) THEN
DO 15 J=MFIT+1,MA
YM=YM-Af(LISTA(J))*AFUNC(LISTA(J))
15 CONTINUE
ENDIF
SIG2I=1./SIG(I)**2
DO 17 J=1,MFIT
WT=AFUNC(LISTA(J))*SIG2I
DO 16 K=1,J
COVAR(J,K)=COVAR(J,K)+WT*AFUNC(LISTA(K))
16 CONTINUE
BETA(J)=BETA(J)+YM*WT
17 CONTINUE
18 CONTINUE
IF (MFIT.GT.1) THEN
DO 21 J=2,MFIT
DO 19 K=1,J-1
COVAR(K,J)=COVAR(J,K)
19 CONTINUE
21 CONTINUE
ENDIF
CALL GAUSSJ(COVAR,MFIT,NCVM,BETA,1,1)
DO 22 J=1,MFIT
Af(LISTA(J))=BETA(J)
22 CONTINUE
CHISQ=0.
DO 24 I=1,NF
CALL FUNCS(xf(I),AFUNC,MA)
SUM=0.
DO 23 J=1,MA
SUM=SUM+Af(J)*AFUNC(J)
23 CONTINUE
CHISQ=CHISQ+((Yf(I)-SUM)/SIG(I))**2
24 CONTINUE
CALL COVSRT(COVAR,NCVM,MA,LISTA,MFIT)
RETURN
END
SUBROUTINE COVSRT(COVAR,NCVM,MA,LISTA,MFIT)
DIMENSION COVAR(NCVM,NCVM),LISTA(MFIT)
DO 12 J=1,MA-1
DO 11 I=J+1,MA
COVAR(I,J)=0.
11 CONTINUE
12 CONTINUE
DO 14 I=1,MFIT-1
DO 13 J=I+1,MFIT
IF(LISTA(J).GT.LISTA(I)) THEN
COVAR(LISTA(J),LISTA(I))=COVAR(I,J)
ELSE
COVAR(LISTA(I),LISTA(J))=COVAR(I,J)
ENDIF
13 CONTINUE
14 CONTINUE
SWAP=COVAR(1,1)
DO 15 J=1,MA
COVAR(1,J)=COVAR(J,J)
COVAR(J,J)=0.
15 CONTINUE
COVAR(LISTA(1),LISTA(1))=SWAP
DO 16 J=2,MFIT
COVAR(LISTA(J),LISTA(J))=COVAR(1,J)
16 CONTINUE
DO 18 J=2,MA
DO 17 I=1,J-1
COVAR(I,J)=COVAR(J,I)
17 CONTINUE
18 CONTINUE
RETURN
END
SUBROUTINE GAUSSJ(A,N,NP,B,M,MP)
PARAMETER (NMAX=100)
DIMENSION A(NP,NP),B(NP,MP),IPIV(NMAX),INDXR(NMAX),INDXC(NMAX)
DO 11 J=1,N
IPIV(J)=0
11 CONTINUE
DO 22 I=1,N
BIG=0.
DO 13 J=1,N
IF(IPIV(J).NE.1)THEN
DO 12 K=1,N
IF (IPIV(K).EQ.0) THEN
IF (ABS(A(J,K)).GE.BIG)THEN
BIG=ABS(A(J,K))
IROW=J
ICOL=K
ENDIF
ELSE IF (IPIV(K).GT.1) THEN
PAUSE '1,Singular matrix'
ENDIF
12 CONTINUE
ENDIF
13 CONTINUE
IPIV(ICOL)=IPIV(ICOL)+1
IF (IROW.NE.ICOL) THEN
DO 14 L=1,N
DUM=A(IROW,L)
A(IROW,L)=A(ICOL,L)
A(ICOL,L)=DUM
14 CONTINUE
DO 15 L=1,M
DUM=B(IROW,L)
B(IROW,L)=B(ICOL,L)
B(ICOL,L)=DUM
15 CONTINUE
ENDIF
INDXR(I)=IROW
INDXC(I)=ICOL
IF (A(ICOL,ICOL).EQ.0.) PAUSE '2,Singular matrix.'
PIVINV=1./A(ICOL,ICOL)
A(ICOL,ICOL)=1.
DO 16 L=1,N
A(ICOL,L)=A(ICOL,L)*PIVINV
16 CONTINUE
DO 17 L=1,M
B(ICOL,L)=B(ICOL,L)*PIVINV
17 CONTINUE
DO 21 LL=1,N
IF(LL.NE.ICOL)THEN
DUM=A(LL,ICOL)
A(LL,ICOL)=0.
DO 18 L=1,N
A(LL,L)=A(LL,L)-A(ICOL,L)*DUM
18 CONTINUE
DO 19 L=1,M
B(LL,L)=B(LL,L)-B(ICOL,L)*DUM
19 CONTINUE
ENDIF
21 CONTINUE
22 CONTINUE
DO 24 L=N,1,-1
IF(INDXR(L).NE.INDXC(L))THEN
DO 23 K=1,N
DUM=A(K,INDXR(L))
A(K,INDXR(L))=A(K,INDXC(L))
A(K,INDXC(L))=DUM
23 CONTINUE
ENDIF
24 CONTINUE
RETURN
END
c
subroutine funcs(ydat,sydat,afun,ma)
c
c
end
Crear cuenta