C @(#)fit_min.for 10.4 (ESO-IPG) 2/9/96 12:10:16 C=========================================================================== C Copyright (C) 1995 European Southern Observatory (ESO) C C This program is free software; you can redistribute it and/or C modify it under the terms of the GNU General Public License as C published by the Free Software Foundation; either version 2 of C the License, or (at your option) any later version. C C This program is distributed in the hope that it will be useful, C but WITHOUT ANY WARRANTY; without even the implied warranty of C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the C GNU General Public License for more details. C C You should have received a copy of the GNU General Public C License along with this program; if not, write to the Free C Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, C MA 02139, USA. C C Corresponding concerning ESO-MIDAS should be addressed as follows: C Internet e-mail: midas@eso.org C Postal address: European Southern Observatory C Data Management Division C Karl-Schwarzschild-Strasse 2 C D 85748 Garching bei Muenchen C GERMANY C=========================================================================== C C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& SUBROUTINE FCN(NFrePm,g,f,x,flag,fuser) C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& c dichiarazioni IMPLICIT DOUBLE PRECISION (a-h,o-z) integer NFrePm ! number of FREE param: ! value may change during minimization integer flag INTEGER NTOT,NTOT2,NPAR,NPAR2 parameter (NTOT=5000,NTOT2=2*NTOT,NPAR=50,NPAR2=100) DOUBLE PRECISION fuser DOUBLE PRECISION gammq,equwid !recipes & routine DOUBLE PRECISION x(NPAR),g(NPAR),f DOUBLE PRECISION tau(NTOT),spe_teo(NTOT),res(NTOT),pixsiz(NTOT) DOUBLE PRECISION alam(NTOT),speoss(NTOT),sig_n(NTOT),stdev(NTOT) DOUBLE PRECISION gaunor(NTOT),fitt(NTOT),chi,prob DOUBLE PRECISION bsum,sum1,sum,chi_rid,chi2,andf integer kk,il,ndf,iok DOUBLE PRECISION ErrPos(NPAR),ErrNeg(NPAR),ErrPar(NPAR) DOUBLE PRECISION Glob(NPAR),w(NPAR2) DOUBLE PRECISION ErLaPo(NPAR2),ErLaNe(NPAR2),ErLaPa(NPAR2) DOUBLE PRECISION ErNnPo(NPAR2),ErNnNe(NPAR2),ErNnPa(NPAR2) DOUBLE PRECISION ErBbPo(NPAR2),ErBbNe(NPAR2),ErBbPa(NPAR2) DOUBLE PRECISION ErBtPo(NPAR2),ErBtNe(NPAR2),ErBtPa(NPAR2) integer ninf(NPAR),nsup(NPAR),k C!!!!!!! C Dati globali (ma limitati a MID_REL_INCL:fit_min) C!!!!!!! DOUBLE PRECISION LamEm(NPAR2),Gam(NPAR2),Fos(NPAR2) DOUBLE PRECISION LamCen(NPAR2),BDopp(NPAR2),Coeff(4,NPAR2) DOUBLE PRECISION BTur(NPAR2),NHFit(NPAR2),FitMn(NPAR2) DOUBLE PRECISION FitMx(NPAR2) integer npunti,NIntv,NRgTot,ij,jj,i integer NTotPm,IPar(4,NPAR2) common/dati/LamEm,Gam,Fos,LamCen,BDopp,BTur,NHFit,FitMn,FitMx, 1 NPunti,NIntv,NRgTot,NTotPm,IPar,Coeff if (flag.eq.1) then C!!!!!!! C Inizializza FCN leggendo la configurazione di righe di default C!!!!!!! c debug c write(6,*)'Start 1 run - read set up' c end debug call InizParam(iok) do i=1,Npunti res(i)=0.D0 end do c debug c write(6,*) LamEm(1),Gam(1),Fos(1),LamCen(1),BDopp(1),NHFit(1) C!!!!!!! C Legge il file con lo spettro, il S/N e la FWHM C!!!!!!! open(10,file='fdummy.spe',status='old') jj=1 10 read(10,*,end=15) alam(jj),pixsiz(jj),speoss(jj),sig_n(jj), + stdev(jj) jj=jj+1 go to 10 15 close(10) npunti=jj-1 call InitConvGauss(stdev,gaunor,npunti) c debug c write(6,*) alam(1), speoss(1), sig_n(1) c write(6,*) alam(2), speoss(2), sig_n(2) c C!!!!!!! C Ricava gli indici relativi agli estremi tra cui calcolare il chi2 C!!!!!!! do ij = 1,NIntv do i=1,NPunti if (FitMn(ij).eq.alam(i)) then ninf(ij)=i else if (FitMn(ij).gt.alam(i).and.FitMn(ij). 1 lt.alam(i+1)) then ninf(ij)=i+1 end if if (FitMx(ij).eq.alam(i)) then nsup(ij)=i else if (FitMx(ij).gt.alam(i).and.FitMx(ij). 1 lt.alam(i+1)) then nsup(ij)=i end if end do end do end if ! fine inizializzazione c write(6,*)alam(1),alam(2),alam(128) C!!!!!!! C OPERAZIONI ORDINARIE C!!!!!!! c debug c write(6,*)'start run' c write(6,*)'lam: ',alam(1),alam(2),alam(128) c write(6,*)'function value: ',f c write(6,*)'Nparam: ',n c write(6,*)'Param: ',x(1),x(2),x(3) c end debug x(NTotPm+1)=0.D0 do i=1,NRGTot LamCen(i)=x(IPar(1,i))*Coeff(1,i) NHFit(i)=x(IPar(2,i))*Coeff(2,i) BDopp(i)=x(IPar(3,i))*Coeff(3,i) BTur(i)=x(IPar(4,i))*Coeff(4,i) end do c write(6,*) LamEm(1),Gam(1),Fos(1),LamCen(1),BDopp(1),NHFit(1), c 1 BTur(1) c do kji =1,NFrePm c write(6,*)x(kji) c end do c write(6,*)' end run' c C!!!!!!! C calcolo della profondita' ottica C!!!!!!! do i=1,NPunti tau(i)=0.D0 end do c debug c write(6,*)NPunti,' Npunti' c write(6,*)'alam:' c write(6,*)alam(1),alam(2),alam(128) c end debug do i=1,NRgTot bsum=dsqrt(BDopp(i)*BDopp(i)+BTur(i)*BTur(i)) call OptDep(LamCen(i),bsum,NHFit(i),LamEm(i), 1 fos(i),gam(i),alam,NPunti,tau) end do c write(6,*)alam(1),alam(2),alam(128) c debug c write(6,*)'tau:' c write(6,*)tau(1),tau(2),tau(128) c write(6,*)' end run' C!!!!!!! C calcolo dello spettro non convoluto C!!!!!!! do i=1,NPunti fitt(i)=dexp(-tau(i)) end do C!!!!!!! C Convolution with Instrumental Response C!!!!!!! call ConvGauss(alam,pixsiz,fitt,stdev,gaunor,spe_teo,NPunti) c do ij=1,NPunti c if (spe_teo(ij).gt.1.D0) spe_teo(ij)=1.D0 c end do c debug c write(6,*)'Spettro teorico:' c write(6,*)spe_teo(1),spe_teo(2),spe_teo(128) c end debug C!!!!!!! c calcolo chi quadro C!!!!!!! kk=0 sum1=0.D0 f=0.D0 do ij = 1,NIntv do i=ninf(ij),nsup(ij) sum=speoss(i)-spe_teo(i) sum1=sum*sum/sig_n(i) f=f+sum1 kk =kk+1 end do end do C!!!!!!! C Fine operazione ORDINARIE (continua solo se ha trovato minimo) C!!!!!!! c debug c write(6,*)'lam: ',alam(1),alam(2),alam(128) c write(6,*)'function value: ',f c write(6,*)'Nparam: ',n c write(6,*)'Param: ',x(1),x(2),x(3) c write(6,*)'End run' c end debug if (flag.ne.3) return c debug c call ConvGauss2(alam,fitt,stdev,gaunor,spe_teo,NPunti) c end debug C!!!!!!! c calcolo probabilita' del fit C!!!!!!! chi=f ndf=kk-NFrePm andf=ndf chi_rid=chi/andf andf=ndf/2.D0 chi2=chi/2.D0 prob=gammq(andf,chi2) c!!!!!!! C Retrieves errors c!!!!!!! do i=1,NTotPm call MNERRS(i,ErrPos(i),ErrNeg(i),ErrPar(i),Glob(i)) end do do i=1,NRGTot ErLaPo(i)=ErrPos(IPar(1,i))*Coeff(1,i) ErLaNe(i)=ErrNeg(IPar(1,i))*Coeff(1,i) ErLaPa(i)=ErrPar(IPar(1,i))*Coeff(1,i) ErNnPo(i)=ErrPos(IPar(2,i))*Coeff(2,i) ErNnNe(i)=ErrNeg(IPar(2,i))*Coeff(2,i) ErNnPa(i)=ErrPar(IPar(2,i))*Coeff(2,i) ErBbPo(i)=ErrPos(IPar(3,i))*Coeff(3,i) ErBbNe(i)=ErrNeg(IPar(3,i))*Coeff(3,i) ErBbPa(i)=ErrPar(IPar(3,i))*Coeff(3,i) ErBtPo(i)=ErrPos(IPar(4,i))*Coeff(4,i) ErBtNe(i)=ErrNeg(IPar(4,i))*Coeff(4,i) ErBtPa(i)=ErrPar(IPar(4,i))*Coeff(4,i) end do C Calcola la larghezza equivalente do il=1,NRgTot bsum=dsqrt(BDopp(il)*BDopp(il)+BTur(il)*BTur(il)) c write(strout,*)LamCen(il),bsum,NHFit(il) c 1 LamEm(il),fos(il),gam(il) c call sttdis(strout,0,il) w(il)=Equwid(LamCen(il),bsum,NHFit(il), 1 LamEm(il),fos(il),gam(il)) end do c!!!!!!! c Aggiorna File fdummy.res c!!!!!!! open(11,file='fdummy.res',status='old',err=33) close(11,status='delete') 33 open(11,file='fdummy.res',status='new') write(11,*)chi_rid,prob do i=1,NRgTot write(11,'(17G19.9E3)')LamCen(i),NHfit(i),BDopp(i),BTur(i), 1 w(i),ErLaPa(i),ErNnPa(i),ErBbPa(i),ErBtPa(i), 2 ErLaPo(i),ErLaNe(i),ErNnPo(i),ErNnNe(i), 3 ErBbPo(i),ErBbNe(i),ErBtPo(i),ErBtNe(i) end do close(11) c!!!!!!! c Aggiorna file fdummy.grf c!!!!!!! open(20,file='fdummy.grf',status='old',err=305) close(20,status='delete') 305 open(20,file='fdummy.grf',status='new') ! calcola i residui do ij = 1,NIntv do i=ninf(ij),nsup(ij) res(i)=speoss(i)-spe_teo(i) end do end do do k=1,npunti c if (spe_teo(k).gt.1.) spe_teo(k)=1.D0 write(20,'(4E16.7)') alam(k),spe_teo(k),res(k),fitt(k) end do close (20) return end FUNCTION H(A,V) C *************** C FUNCTION DI VOIGT IMPLICIT DOUBLE PRECISION (A-H,O-Z) LOGICAL Q VV=V*V Q=A.LT.0.2 IF(Q.AND.V.GT.5.)GOTO 1 IF(.NOT.Q.AND.(A.GT.1.4.OR.A+V.GT.3.2))GOTO 2 HO=DEXP(-VV) H2=(1.-2.*VV)*HO IF(V.GT.2.4)GOTO 3 IF(V.GT.1.3)GOTO 4 H1=(.42139D0*VV-2.34358D0*V+3.28868D0)*VV-.15517D0*V-1.1247D0 5 IF(Q)GOTO 6 HH1=H1+HO*1.12838D0 HH2=H2+HH1*1.12838D0-HO HH3=(1.-H2)*.37613D0-HH1*.66667D0*VV+HH2*1.12838D0 HH4=(3.D0*HH3-HH1)*.37613D0+HO*.66667D0*VV*VV H=((((HH4*A+HH3)*A+HH2)*A+HH1)*A+HO)*(((-.122727278D0*A 1 +.532770573D0)*A-.96284325D0)*A+.979895032D0) RETURN 1 H=((2.12D0/VV+.8463D0)/VV+.5642D0)*A/VV RETURN 2 AA=A*A U=(AA+VV)*1.4142D0 UU=U*U H=((((AA-10.D0*VV)*AA*3.D0+15.D0*VV*VV)/UU+3.D0*VV-AA)/UU+1.D0)* 1 A*.79788D0/U RETURN 3 H1=((-.0032783D0*VV+.0429913D0*V-.188326D0)*VV+.278712D0* > V+.55415D0)/(VV-1.5D0) GOTO 5 4 H1=(-.220416D0*VV+1.989196D0*V-6.61487D0)*VV+9.39456D0*V-4.4848D0 GOTO 5 6 H=(H2*A+H1)*A+HO RETURN END C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& function Equwid(LC,B,N,LE,fos,gam) C C Calcola la larghezza equivalente di una riga C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION LC,B,LE,N,fos,gam DOUBLE PRECISION le0 integer nint double precision offset,tau common /datiriga/al,bd,ah,le0,fos0,gam0 external GG le0=le fos0=fos gam0=gam al=LC bd=B ah=N offset = 4. 10 offset = offset +2. alambda2 = LC+offset tau=0. call OptDep(LC,B,N,LE,fos,gam,alambda2,1,tau) if (tau.gt.1.D-3) goto 10 alambda2=LC+offset alambda1=LC-offset nint=2.*offset/0.05 +1 call sim(alambda1,alambda2,sint,nint,GG) Equwid=sint/LC*LE return end C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& FUNCTION GG(ALAMM) C C Funzione ausiliaria per il calcolo della larghezza equiv. C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION LC,B,LE,N,fos,gam INTEGER NPUNTI common /datiriga/lc,b,n,le,fos,gam tau=0. NPunti=1 call OptDep(LC,B,N,LE,fos,gam,alamm,Npunti,tau) GG=1.D0-dexp(-tau) return end C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& SUBROUTINE SIM(AA,BB,SOMMA,NINT,G) C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER NINT,JK,NN EXTERNAL G BA=BB-AA SOMMA=G(AA)+G(BB) VINT=BA/NINT XD=AA+VINT/2.D0 SD=G(XD) SP=0. NN=NINT-1 DO 1 JK=1,NN XP=XD+VINT/2.D0 XD=XP+VINT/2.D0 SP=SP+G(XP) 1 SD=SD+G(XD) SOMMA=SOMMA+2.D0*SP+4.D0*SD SOMMA=SOMMA*VINT/6.D0 RETURN END C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& subroutine InizParam(flag) C C Legge una configurazione di righe dal file fitly.fcn C e genera i vettori di parametri delle righe C che usa FCN C Ritorna una flag di errore C C flag: integer flag di errore: se 0 ok, C se <0 errore C C C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& integer flag C IMPLICIT none INTEGER NTOT,NTOT2,NPAR,NPAR2 parameter (NTOT=5000,NTOT2=NTOT*2,NPAR=50,NPAR2=100) DOUBLE PRECISION LamEm(NPAR2),Gam(NPAR2),Fos(NPAR2) DOUBLE PRECISION LamCen(NPAR2),BDopp(NPAR2),Coeff(4,NPAR2) DOUBLE PRECISION BTur(NPAR2),NHFit(NPAR2),FitMn(NPAR2) DOUBLE PRECISION FitMx(NPAR2) integer npunti,NIntv,NRgTot integer NTotPm,IPar(4,NPAR2) common/dati/LamEm,Gam,Fos,LamCen,BDopp,BTur,NHFit,FitMn,FitMx, 1 NPunti,NIntv,NRgTot,NTotPm,IPar,Coeff C character*20 Junk integer i,k,j C&&&&&&& C Azzera un po' di variabili.. C&&&&&&& NTotPm=0 do i=1,NPAR2 LamEm(i)=0.D0 Gam(i)=0.D0 Fos(i)=0.D0 LamCen(i)=0.D0 BDopp(i)=0.D0 BTur(i)=0.D0 NHFit(i)=0.D0 FitMn(i)=0.D0 FitMx(i)=0.D0 do k=1,4 Coeff(k,i)=0.D0 IPar(k,i)=0 end do end do flag=-1 C&&&&&&& C& Apre file C&&&&&&& open(10,file='fdummy.fcn',iostat=i,status='old') if (i.ne.0) then flag=-1 close(10) return end if C&&&&&&& C& Legge parametri innesco C&&&&&&& 101 format(4I3,G16.8,F4.1,G16.8,F4.1,F11.4,F9.5,G13.5) read(10,*)NRgTot do j=1,NRgTot 90 read(10,1001,err=99)(IPar(i,j),i=1,4),(Coeff(i,j),i=1,4) 1 ,LamEm(j),Fos(j),Gam(j) end do 1001 format(4I3,7G16.8) C gets number of parameters do j=1,NRgTot do i=1,4 if (IPar(i,j).gt.NTotPm) NTotPm= IPar(i,j) end do end do C&&&&&&& C Reads limits C&&&&&&& read(10,*,end=99,err=99)NIntv do i = 1,NIntv read(10,*,end=99,err=99)FitMn(i),FitMx(i) end do flag=0 99 continue close(10) return end C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& subroutine Minmze C C Chiama la routine di minimizzazione MINUIT C C C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& external FCN integer ista,iwrite,iread,ipunch c include 'MID_REL_INCL:fit_var.inc' open(90,file='fdummy.min',status='old',err=93) open(91,file='punch.dat',status='old',err=10) close(91,status='delete') 10 open(91,file='punch.dat',status='new',iostat=ista) open(92,file='fdummy.jou',status='old',err=20) close(92,status='delete') 20 open(92,file='fdummy.jou',status='new',iostat=ista) iwrite=92 iread=90 ipunch=91 call mintio(iread,iwrite,ipunch) CALL MINUIT(FCN,0) close(90) close(91,status='delete') close(92) 93 return end C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& Subroutine InitConvGauss(stdev,gaunor,npunti) C C Prepare variables for ConvGauss from fwhm C C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& INTEGER NTOT,NTOT2,NPAR,NPAR2 parameter (NTOT=5000,NTOT2=2*NTOT,NPAR=50,NPAR2=100) integer npunti,ijk DOUBLE PRECISION stdev(npunti),gaunor(npunti) DOUBLE PRECISION pi pi=dacos(-0.1D+01) c DOUBLE PRECISION ax,step c DOUBLE PRECISION gausstab(NTOT2) c common/gauss/gausstab,step,nstep c nstep=NTOT c step = 3.5/nstep c do ij=1,nstep c rij=ij - 1 c ax = rij * step c gausstab(ij)= dexp(-ax*ax) c end do DO IJK=1,npunti stdev(ijk) = stdev(ijk)/2.D0/dsqrt(2.D0*dlog(2.D0)) ! FWHM -> sigma gaunor(ijk) = 1.D0 / dsqrt(2.D0 * pi) /stdev(ijk) ! gauss. norm. stdev(ijk) = stdev(ijk) * dsqrt(2.D0) ! sigma-> s*rad2 end do return end C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& Subroutine ConvGauss(alam,pixsiz,fitt,sig,norm,spe_t,NPunti) C C Performs convolution of a spectrum with a gaussian response c function variyng with wavelength c c alam: spectrum wavelengths c fitt: spectrum data c sig: sqrt(2.) * stdev gaussian (stdev = FWHM / 2/sqrt(2*alog(2)) ) c norm: stdev * sqrt(2*3.14) (gaussian normalization) c spe_t: convolved spectrum C C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& C IMPLICIT NONE integer NPunti DOUBLE PRECISION alam(Npunti),fitt(NPunti),sig(NPunti) DOUBLE PRECISION spe_t(NPunti),norm(NPunti),pixsiz(NPunti) DOUBLE PRECISION gauss c integer nstep,i,ind c DOUBLE PRECISION gausstab(10000),step c common/gauss/gausstab,step,nstep integer il,ic DOUBLE PRECISION sum,dlam do il=1,NPunti sum=0.D0 do ic=il,Npunti dlam = (alam(ic) - alam(il)) / sig(ic) c ind = idnint(dabs(dlam)/step) + 1 c if (ind.gt.nstep) goto 100 c sum = sum + norm(ic)*gausstab(ind)*fitt(ic) gauss = dexp(-dlam * dlam) if (gauss.lt.1D-5) goto 100 sum = sum + norm(ic)*gauss*fitt(ic) end do 100 continue if (gauss.gt.1D-5) then dlam = dlam+pixsiz(il)/ sig(Npunti) gauss = dexp(-dlam * dlam) sum = sum + norm(Npunti)*gauss*1.D0 !convolvo continuo goto 100 endif do ic=il-1,1,-1 dlam = (alam(ic) - alam(il))/sig(ic) c ind = idnint(dabs(dlam)/step) + 1 c if (ind.gt.nstep) goto 200 c sum = sum + norm(ic)*gausstab(ind)*fitt(ic) gauss = dexp(-dlam * dlam) if (gauss.lt.1D-5) goto 200 sum = sum + norm(ic)*gauss*fitt(ic) end do 200 continue if (il.eq.1) then gauss =1.D0 dlam = 0.D0 endif 201 if (gauss.gt.1D-5) then dlam = dlam-pixsiz(il)/ sig(1) gauss = dexp(-dlam * dlam) sum = sum + norm(1)*gauss*1.D0 !convolvo continuo goto 201 endif ! moltiplica per il pixel size. spe_t(il) = sum * pixsiz(il) end do return end C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& SUBROUTINE OptDep(cla,b,nh,lamem,f,gamma,lam,nlam,tau) C C Calcola la profondita' ottica relativa ad una riga e la SOMMA C a quella gia' esistente. C Versione FAST: il calcolo e' limitato alla zona in cui tau>10^-7 C C Variabili: C cla DOUBLE PRECISION Input lambda osservata della riga C b DOUBLE PRECISION " b Doppler riga C nh DOUBLE PRECISION " NH riga C lamem DOUBLE PRECISION " lamda emissione C f DOUBLE PRECISION " Forza oscillatore C gamma DOUBLE PRECISION " Damping constant C lam DOUBLE PRECISION(nlam) " vettore lambda C nlam INTEGER " numero punti C tau DOUBLE PRECISION(nlam) In/Out vettore prof. ottica C C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION cla,b,nh,f,lamem,gamma,lam(nlam),tau(nlam) DOUBLE PRECISION SOGLIA PARAMETER (SOGLIA=1.D-7) INTEGER NLAM,I,ind0 pi=dacos(-0.1D+01) sp=dsqrt(pi) r0=2.81794D-13 c=3.0D+05 a=0.D0 hh=0.D0 sig1=0.D0 sig2=0.D0 dlamD=b*cla/c a=gamma*lamem/(4.D0*pi*b*1.D13) sig1=r0*lamem*1.d-8*f*sp/dlamD*cla do i=1,nlam if (lam(i).gt.cla) goto 10 end do 10 ind0=i do i=ind0,nlam v=dabs((cla-lam(i))/dlamD) hh=h(a,v) sig2=sig1*hh tau(i)=10.D0**nh*sig2+tau(i) if (tau(i).lt.SOGLIA) goto 20 end do 20 continue do i=ind0-1,1,-1 v=dabs((cla-lam(i))/dlamD) hh=h(a,v) sig2=sig1*hh tau(i)=10.D0**nh*sig2+tau(i) if (tau(i).lt.SOGLIA) goto 30 end do c debug write(6,*)lam(1),lam(nlam) 30 continue return end C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& subroutine CHECK(gptau,NomEle,lambda,b,NH,fullwh,nrighe,x, 1 pixs,npt,flux,unconv) C C Compute an absorption spectrum C C C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& implicit none include 'MID_REL_INCL:fit_var.inc' ! for atompar.dat integer npt,nrighe character*14 NomEle(nrighe) double precision gptau,lambda(nrighe),b(nrighe) double precision NH(nrighe),fullwh(nrighe) double precision x(npt),flux(npt),unconv(npt),pixs(npt) double precision lamem,f,gamma double precision gaunor(NMXSPE),tau(NMXSPE) double precision stdev(NMXSPE),gptau0 integer j,i c integer pp do i=1,npt stdev(i)=fullwh(i) end do gptau0=dexp(-gptau) call InitConvGauss(stdev,gaunor,npt) C&&&&&&& C Calcola il profilo teorico C&&&&&&& do i=1,npt tau(i)=0.0 end do do i=1,nrighe do j = 1,AT_N if (NomEle(i).eq.AT_NAM(j)) then lamem = AT_LAM(j) f = AT_FOS(j) gamma = AT_GAM(j) goto 456 end if end do 456 continue call OptDep(lambda(i),b(i),NH(i),lamem, 1 f,gamma,x,npt,tau) end do do i=1,npt unconv(i)=exp(-tau(i)) * gptau0 end do call ConvGauss(x,pixs,unconv,stdev,gaunor,flux,Npt) c debug c do i=1,npt c if (flux(i).lt.-0.5.or.flux(i).gt.1.5) then c open(20,file='test2.dbx',status='unknown') c write(20,*)'lambda, flusso:, flusso conv.' c do pp = 1,npt c write(20,*)x(pp),tau(pp),flux(pp) c end do c close(20) c stop c end if c end do c end debug return end C****** function INTRAC() C C required by MINUIT C****** logical INTRAC INTRAC=.FALSE. return end CC C This is a small subset of the "Numerical Recipes" routines, C currently used by the program. They are exactly identical to C those of the original library, except for a different warning C (instead of 'A too large, ITMAX too small' , a C Midas warning is inserted) FUNCTION GAMMQ(A,X) IMPLICIT DOUBLE PRECISION(A-H,O-Z) IMPLICIT INTEGER (I-N) IF(X.LT.0..OR.A.LE.0.)PAUSE IF(X.LT.A+1.)THEN CALL GSER(GAMSER,A,X,GLN) GAMMQ=1.D0-GAMSER ELSE CALL GCF(GAMMCF,A,X,GLN) GAMMQ=GAMMCF ENDIF RETURN END SUBROUTINE GSER(GAMSER,A,X,GLN) IMPLICIT INTEGER (I-N) IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER (ITMAX=100,EPS=3.E-7) GLN=GAMMLN(A) IF(X.LE.0.)THEN IF(X.LT.0.)PAUSE GAMSER=0. RETURN ENDIF AP=A SUM=1./A DEL=SUM DO 11 N=1,ITMAX AP=AP+1. DEL=DEL*X/AP SUM=SUM+DEL IF(DABS(DEL).LT.DABS(SUM)*EPS)GO TO 1 11 CONTINUE c PAUSE 'A too large, ITMAX too small' call WrnMsg('Prob(chi2) may be inaccurate') 1 GAMSER=SUM*DEXP(-X+A*DLOG(X)-GLN) RETURN END SUBROUTINE GCF(GAMMCF,A,X,GLN) IMPLICIT INTEGER (I-N) IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER (ITMAX=100,EPS=3.E-7) GLN=GAMMLN(A) GOLD=0. A0=1. A1=X B0=0. B1=1. FAC=1. DO 11 N=1,ITMAX AN=DFLOAT(N) ANA=AN-A A0=(A1+A0*ANA)*FAC B0=(B1+B0*ANA)*FAC ANF=AN*FAC A1=X*A0+ANF*A1 B1=X*B0+ANF*B1 IF(A1.NE.0.D0)THEN FAC=1.D0/A1 G=B1*FAC IF(DABS((G-GOLD)/G).LT.EPS)GO TO 1 GOLD=G ENDIF 11 CONTINUE c PAUSE 'A too large, ITMAX too small' call WrnMsg('Prob(chi2) may be inaccurate') 1 GAMMCF=DEXP(-X+A*DLOG(X)-GLN)*G RETURN END SUBROUTINE PIKSR2(N,ARR,BRR) IMPLICIT DOUBLE PRECISION(A-H,O-Z) IMPLICIT INTEGER (I-N) DIMENSION ARR(N),BRR(N) DO 12 J=2,N A=ARR(J) B=BRR(J) DO 11 I=J-1,1,-1 IF(ARR(I).LE.A)GO TO 10 ARR(I+1)=ARR(I) BRR(I+1)=BRR(I) 11 CONTINUE I=0 10 ARR(I+1)=A BRR(I+1)=B 12 CONTINUE RETURN END FUNCTION GAMMLN(XX) IMPLICIT DOUBLE PRECISION(A-H,O-Z) IMPLICIT INTEGER (I-N) DOUBLE PRECISION COF(6),STP,HALF,ONE,FPF,X,TMP,SER DATA COF,STP/76.18009173D0,-86.50532033D0,24.01409822D0, * -1.231739516D0,.120858003D-2,-.536382D-5,2.50662827465D0/ DATA HALF,ONE,FPF/0.5D0,1.0D0,5.5D0/ X=XX-ONE TMP=X+FPF TMP=(X+HALF)*DLOG(TMP)-TMP SER=ONE DO 11 J=1,6 X=X+ONE SER=SER+COF(J)/X 11 CONTINUE GAMMLN=TMP+DLOG(STP*SER) RETURN END function gasdev(idum) double precision gasdev,v1,v2,r,fac,gset,ran1n external ran1n integer iset,idum iset=0 if(iset.eq.0)then 1 v1=2.D0*ran1n(idum)-1.D0 v2=2.D0*ran1n(idum)-1.D0 r=v1**2.D0+v2**2.D0 if(r.ge.1.)go to 1 fac=dsqrt(-2.*dlog(r)/r) gset=v1*fac gasdev=v2*fac iset=1 else gasdev=gset iset=0 endif return end FUNCTION RAN1N(IDUM) double PRECISION RAN1N DOUBLE PRECISION R(97) INTEGER M1,M2,M3,IA1,IA2,IA3,IC1,IC2,IC3,IFF INTEGER IDUM,ix1,ix2,ix3,j DOUBLE PRECISION RM1,RM2 PARAMETER (M1=259200,IA1=7141,IC1=54773,RM1=3.8580247D-6) PARAMETER (M2=134456,IA2=8121,IC2=28411,RM2=7.4373773D-6) PARAMETER (M3=243000,IA3=4561,IC3=51349) data IFF /1/ IF (IDUM.LT.0.OR.IFF.EQ.0) THEN IFF=1 IX1=MOD(IC1-IDUM,M1) IX1=MOD(IA1*IX1+IC1,M1) IX2=MOD(IX1,M2) IX1=MOD(IA1*IX1+IC1,M1) IX3=MOD(IX1,M3) DO 11 J=1,97 IX1=MOD(IA1*IX1+IC1,M1) IX2=MOD(IA2*IX2+IC2,M2) R(J)=(DFLOAT(IX1)+DFLOAT(IX2)*RM2)*RM1 11 CONTINUE IDUM=1 ENDIF IX1=MOD(IA1*IX1+IC1,M1) IX2=MOD(IA2*IX2+IC2,M2) IX3=MOD(IA3*IX3+IC3,M3) J=1+(97*IX3)/M3 IF(J.GT.97.OR.J.LT.1)PAUSE RAN1N=R(J) R(J)=(DFLOAT(IX1)+DFLOAT(IX2)*RM2)*RM1 RETURN END