C @(#)dblend.for 17.1.1.1 (ES0-DMD) 01/25/02 17:56:13 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 program DBLEND C C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C program DBLEND C C Uni-Sternwarte Goettingen C Adaptation of TVR program TVDBLE, DBLEND command C C purpose: C C This routine provides de-blending of spectral lines by means of a C simultaneous fit of up to six profiles to a selected spectral C region. C C Run by procedure DBLEND C C------------------------------------------------------------------------------ implicit none C------------------------------------------------------------------------------ C C Variable declarations C integer MADRID !mapped section pointer function C INTEGER ACOL !number of words allocated per table row INTEGER ACTVAL !actual number of values returned INTEGER AROW !number of rows allocated in table REAL CENT(6) !centers of lines REAL CMAX(6) !central line intensity INTEGER COLS(3) !numbers of columns to be read from table REAL CONTIN(4) !x,y of 2 continuum points chosen by user character*64 CUNIT !units REAL DATA(50000)!spectrum to be worked on REAL FITLIM(2) !limits of area to be fit character*8 FNAME !file names for individual line fits character*3 FOPT !fit option character*1 FOPT1,FOPT2,FOPT3 !individual characters of FOPT REAL FWHM !full width half max of line character*72 IDENT !ASCII identification character*72 IDNT !Temporary buffer character*60 INFILE !name of input spectrum to be fit character*60 INTAB !name of input table of line parameters INTEGER INTID !table identifier assigned to input table file INTEGER KNUL !no. of null values INTEGER KUN !units character*1 LNO !character representation of line number REAL LSIDE !left FWHM point of line INTEGER LTAB !number of columns in output table INTEGER MODE,MODE1 !fit option numerical code REAL MPAR(18) !line parameters read from table INTEGER NAXIS !dimensionality of image INTEGER NCOL !number of columns in table INTEGER NDATA !number of points in spectrum INTEGER NLINES !number of profiles to be fit INTEGER NO !file number assigned to frame INTEGER NPIX(2) !length of spectrum INTEGER NROW !number of rows in table INTEGER NSC !sorted column in table logical*1 NULL(3) !null flag character*60 OUTFILE !name of output file where fit is to be stored character*60 OUTTAB !name of output table of line parameters INTEGER OUTTID !table identifier assigned to output table file REAL PAR(10) !line parameters to be written to table INTEGER*8 PNTR !pointer to mapped data area of input image INTEGER RCOLS(10) !numbers of columns to be written to table REAL RSIDE !right FWHM point of line REAL SCENT(6) !estimated error of fit center REAL SIGMA(6) !Gaussian sigma of line REAL SMAX(6) !estimated error of fit maximum REAL SPAR(18) !line parameters read from table REAL SSIGMA(6) !estimated error of fit sigma DOUBLE PRECISION START(2) !start coordinate INTEGER STAT !return status for subroutine calls DOUBLE PRECISION STEP(2) !step size character*80 STRING !holder for output messages C INTEGER II,JJ !dummy counter INTEGER NN !dummy value holder C C------------------------------------------------------------------------------ common /VMR/ MADRID(1) C------------------------------------------------------------------------------ data COLS/1,3,5/ data RCOLS/1,2,3,4,5,6,7,8,9,10/ C------------------------------------------------------------------------------ C C Start up program as part of MIDAS environment C CALL stspro('DBLEND') C C Get input parameters C CALL stkrdc('INFRAM',1,1,60,ACTVAL,INFILE,KUN,KNUL,STAT) CALL stkrdc('IN_B',1,1,60,ACTVAL,INTAB,KUN,KNUL,STAT) CALL stkrdc('OUT_A',1,1,60,ACTVAL,OUTFILE,KUN,KNUL,STAT) CALL stkrdc('OUT_B',1,1,60,ACTVAL,OUTTAB,KUN,KNUL,STAT) CALL stkrdc('FITOPT',1,1,3,ACTVAL,FOPT,KUN,KNUL,STAT) CALL upcas(FOPT,FOPT) FOPT1 = FOPT(1:1) FOPT2 = FOPT(2:2) FOPT3 = FOPT(3:3) C C translate fit option to numerical code, with all parameters free as default C if (FOPT1.eq.'W') then MODE = 2 else if (FOPT1.eq.'F') then MODE = 3 else if (FOPT1.eq.'S') then MODE = 4 else MODE = 1 end if C C More weighting options: add a paragraph here C if (FOPT2.eq.'P') then STRING = ' Weighting option chosen' CALL sttput(STRING,STAT) MODE = MODE + 10 end if C C More curve options: add a paragraph here C if (FOPT3.eq.'L') then STRING = ' Lorentz curve chosen' CALL sttput(STRING,STAT) MODE = MODE + 100 end if C C Open input image C CALL STIGET(INFILE,10,0,1,2,NAXIS,NPIX,START,STEP,IDENT, & CUNIT,PNTR,NO,STAT) C if (NAXIS.eq.1) then NDATA = NPIX(1) CALL read1d(MADRID(PNTR),NPIX(1),DATA) else if (NAXIS.eq.2.and.(NPIX(1).eq.1.or.NPIX(2).eq.1)) then if (NPIX(1).eq.1) then NDATA = NPIX(2) START(1) = START(2) STEP(1) = STEP(2) else NDATA = NPIX(1) end if CALL read2d(MADRID(PNTR),NPIX,DATA,NDATA) else STRING = 'Input image is not a suitable type' CALL STTPUT(STRING,STAT) go to 9999 end if C C Open input table C CALL TBTOPN(INTAB,2,INTID,STAT) CALL TBIGET(INTID,NCOL,NROW,NSC,ACOL,AROW,STAT) if (NROW.gt.22) then STRING = 'WARNING: only the first 6 lines in ' // + 'the table will be used' CALL STTPUT(STRING,STAT) end if NN = (NROW-4) / 3 NLINES = min(NN,6) C C Read parameters from table C CALL TBERDR(INTID,1,3,FITLIM(1),NULL,STAT) CALL TBERDR(INTID,2,3,FITLIM(2),NULL,STAT) CALL TBERDR(INTID,3,3,CONTIN(1),NULL,STAT) CALL TBERDR(INTID,4,3,CONTIN(2),NULL,STAT) CALL TBERDR(INTID,3,2,CONTIN(3),NULL,STAT) CALL TBERDR(INTID,4,2,CONTIN(4),NULL,STAT) do II = 5,NROW CALL TBERDR(INTID,II,1,SPAR(II-4),NULL,STAT) CALL TBERDR(INTID,II,2,MPAR(II-4),NULL,STAT) end do JJ = 0 do II = 1,NLINES JJ = JJ + 1 LSIDE = SPAR(JJ) JJ = JJ + 1 CENT(II) = ((SPAR(JJ) - real(START(1))) / real(STEP(1)))+1. CMAX(II) = MPAR(JJ) JJ = JJ + 1 RSIDE = SPAR(JJ) FWHM = abs(RSIDE - LSIDE) / real(STEP(1)) if (MODE.ge.100) then SIGMA(II) = FWHM !Lorentz profiles LTAB = 8 else SIGMA(II) = FWHM / 2.354 !Gaussian profiles LTAB = 10 end if if (JJ.gt.NROW) go to 9999 end do C CALL TBTCLO(INTID,STAT) C C CALL the subroutine that prepares the data points and calls the fit routine C MODE1 = MODE CALL prepfit(DATA,NDATA,FITLIM,CONTIN,NLINES,MODE, & CENT,SIGMA,CMAX,SCENT,SSIGMA,SMAX) MODE = MODE1 C C Create output table C CALL TBTINI(OUTTAB,1,1,LTAB,NLINES,OUTTID,STAT) C C Initialize the columns C CALL TBCINI(OUTTID,10,1,'E12.4',CUNIT,'CENTER',NN,STAT) CALL TBCINI(OUTTID,10,1,'E12.4',CUNIT,'ERR_CENT',NN,STAT) CALL TBCINI(OUTTID,10,1,'E12.4','FLUX','MAXIMUM',NN,STAT) CALL TBCINI(OUTTID,10,1,'E12.4','FLUX','ERR_MAX',NN,STAT) if (LTAB.eq.10) then CALL TBCINI(OUTTID,10,1,'E12.4',CUNIT,'SIGMA',NN,STAT) CALL TBCINI(OUTTID,10,1,'E12.4',CUNIT,'ERR_SIG',NN,STAT) end if CALL TBCINI(OUTTID,10,1,'E12.4',CUNIT,'FWHM',NN,STAT) CALL TBCINI(OUTTID,10,1,'E12.4',CUNIT,'ERR_FWHM',NN,STAT) CALL TBCINI(OUTTID,10,1,'E12.4','FLUX','INT_FLUX',NN,STAT) CALL TBCINI(OUTTID,10,1,'E12.4','FLUX','ERR_SIG',NN,STAT) C C Write resulting fit parameters to output table C do II = 1,NLINES PAR(1) = (CENT(II)-1)*real(STEP(1)) + real(START(1)) PAR(2) = SCENT(II)*real(STEP(1)) PAR(3) = CMAX(II) PAR(4) = SMAX(II) PAR(5) = SIGMA(II)*real(STEP(1)) PAR(6) = SSIGMA(II)*real(STEP(1)) if (LTAB.eq.10) then PAR(7) = PAR(5) * 2.354 PAR(8) = PAR(6) * 2.354 PAR(9) = SQRT(2*3.14159) * PAR(5) * PAR(3) PAR(10) = PAR(9) * ((PAR(4)/PAR(3)) + (PAR(6)/PAR(5))) else PAR(7) = (3.14159/2) * PAR(5) * PAR(3) PAR(8) = PAR(7) * ((PAR(4)/PAR(3))+(PAR(6)/PAR(5))) end if CALL tbrwrr(OUTTID,II,LTAB,RCOLS,PAR,STAT) end do C CALL TBTCLO(OUTTID,STAT) C C Finished with input file C CALL STFCLO(NO,STAT) C C Now to generate output files C NPIX(1) = abs(FITLIM(2) - FITLIM(1)) + 1 START(1) = START(1) + (dble(FITLIM(1)-1) * STEP(1)) IDENT = 'DBLEND fit of spectrum '//INFILE STRING = 'generating file '//OUTFILE CALL sttput(STRING,STAT) CALL STIPUT(OUTFILE,10,1,1,1,NPIX(1),START(1),STEP(1),IDENT, & CUNIT,PNTR,NO,STAT) CALL fitgen(MADRID(PNTR),NPIX(1),FITLIM,CONTIN,NLINES, & CENT,SIGMA,CMAX,MODE) CALL stfclo(NO,STAT) C C Individual lines C if (NLINES.gt.1) then do II = 1,NLINES write (LNO,fmt='(i1)') II IDNT = IDENT//', line number '//LNO IDENT = IDNT FNAME = OUTFILE(1:3)//'L000'//LNO STRING = 'generating file '//FNAME CALL sttput(STRING,STAT) CALL stiput(FNAME,10,1,1,1,NPIX(1),START(1),STEP(1),IDENT, & CUNIT,PNTR,NO,STAT) CALL fitgen(MADRID(PNTR),NPIX(1),FITLIM,CONTIN,1,CENT(II), & SIGMA(II),CMAX(II),MODE) CALL stfclo(NO,STAT) end do end if C C All finished! C 9999 CALL stsepi end subroutine prepfit(IN,NPIX,FITLIM,CONTIN,NLINES,MODE, & CENT,SIGMA,CMAX,SCENT,SSIGMA,SMAX) implicit none INTEGER NLINES !number of profiles to be fit INTEGER NPIX !size of image DOUBLE PRECISION A(40) !parameters to be passed to fitting routine DOUBLE PRECISION A0(40) !iteration limits for parameters A REAL CENT(NLINES) !centers of lines REAL CMAX(NLINES) !central line intensity REAL CONTIN(4) !x,y of 2 continuum points chosen by user DOUBLE PRECISION D(40) !step size for calculating derivatives REAL FITLIM(2) !limits of area to be fit REAL IN(NPIX) !input image INTEGER LIM0 !zeroth pixel in subsection to be fit INTEGER MODE !fit mode INTEGER NPIXS !number of pixels in subsection to be fit DOUBLE PRECISION S ! ?? cf snlfit REAL SCENT(NLINES) !estimated error of fit center DOUBLE PRECISION SIG(40) !error estimates of parameters A REAL SIGMA(NLINES) !Gaussian sigma of line DOUBLE PRECISION SLOPE !slope of continuum REAL SMAX(NLINES) !estimated error of fit maximum REAL SSIGMA(NLINES) !estimated error of fit sigma INTEGER STAT !return status indicator character*80 STRING !character screen for terminal output INTEGER TPIX !pixel counter DOUBLE PRECISION VAR ! ?? cf snlfit DOUBLE PRECISION WAV !average weight of data points DOUBLE PRECISION WW(50000) !weights for data points DOUBLE PRECISION XCONT(2),YCONT(2) !x,y for chosen continuum points DOUBLE PRECISION XX(50000),YY(50000) !data values for least squares fit DOUBLE PRECISION YBACK !value of continuum at given pixel C DOUBLE PRECISION DX,DY !temporary value holders INTEGER II !dummy counter C external GAUSS,LORENTZ data WW /50000*1./ C C Read in portion of spectrum to be fit C NPIXS = nint(abs(FITLIM(2) - FITLIM(1))) + 1 LIM0 = nint(FITLIM(1)) - 1 do II = 1,NPIXS TPIX = LIM0 + II XX(II) = dble(TPIX) YY(II) = dble(IN(TPIX)) if (YY(II).eq.0) then STRING = ' zero data value at pixel ' write (STRING(27:80),fmt='(i4)') II CALL sttput(STRING,STAT) end if end do C C More weighting options: add a paragraph here C if (MODE.ge.10) then do II = 1,NPIXS WW(II) = 1./abs(YY(II)) if (WW(II).eq.0) then STRING = ' zero weight at pixel ' write (STRING(23:80),fmt='(i4)') II CALL sttput(STRING,STAT) end if end do end if C C Subtract continuum C do II = 1,2 XCONT(II) = dble(CONTIN(II)) YCONT(II) = dble(CONTIN(II+2)) end do SLOPE = (YCONT(2) - YCONT(1)) / (XCONT(2) - XCONT(1)) do II = 1,NPIXS DX = XX(II) - XCONT(1) DY = DX * SLOPE YBACK = YCONT(1) + DY YY(II) = YY(II) - YBACK if (YY(II).eq.0) then STRING = ' zero data value at pixel ' write (STRING(27:80),fmt='(i4)') II CALL sttput(STRING,STAT) end if end do C C Fill variables in common block C do II = 1,40 A(II) = 0. A0(II) = 0. SIG(II) = 0. D(II) = 1.0d-6 end do do II = 1,NLINES A(II) = dble(CENT(II)) A(II+NLINES) = dble(SIGMA(II)) DX = dble(CENT(II)) - XCONT(1) DY = DX * SLOPE YBACK = YCONT(1) + DY A(II+2*NLINES) = dble(CMAX(II)) - YBACK A0(II) = 0.01 * dble(SIGMA(II)) A0(II+NLINES) = A0(II) A0(II+2*NLINES) = 0.01 * dble(CMAX(II)) end do if (MODE.ge.100) then CALL snlfit(MODE,XX,YY,WW,NPIXS,NLINES,LORENTZ, & A,A0,SIG,D,WAV,VAR,S) else CALL snlfit(MODE,XX,YY,WW,NPIXS,NLINES,GAUSS, & A,A0,SIG,D,WAV,VAR,S) end if C C Re-parse variables in common block C do II = 1,NLINES CENT(II) = real(A(II)) SIGMA(II) = real(A(II+NLINES)) CMAX(II) = real(A(II+2*NLINES)) SCENT(II) = real(SIG(II)) SSIGMA(II) = real(SIG(II+NLINES)) SMAX(II) = real(SIG(II+2*NLINES)) end do return end subroutine fitgen(OUT,NPIX,FITLIM,CONTIN,NLINES,CENT,SIGMA, & CMAX,MODE) implicit none INTEGER NLINES !number of profiles to be fit INTEGER NPIX !size of image DOUBLE PRECISION AA(40) !parameter array for fit subroutines REAL CENT(NLINES) !centers of lines REAL CMAX(NLINES) !central line intensity REAL CONTIN(4) !x,y of 2 continuum points chosen by user REAL FITLIM(2) !limits of area to be fit REAL GSUM !calculated continuum background at a given pixel INTEGER MODE !code for type of fit to be done REAL OUT(NPIX) !output image INTEGER PIX REAL SIGMA(NLINES) !Gaussian sigma of line REAL SLOPE !slope of continuum DOUBLE PRECISION X0,Y0 !a point on the fitted curve REAL XCONT(2),YCONT(2) !x,y for chosen continuum points REAL DX,DY !temporary value holders INTEGER II,JJ !dummy counter do II = 1,2 XCONT(II) = CONTIN(II) YCONT(II) = CONTIN(II+2) end do SLOPE = (YCONT(2) - YCONT(1)) / (XCONT(2) - XCONT(1)) C STRING = ' FITLIM = ' C write (STRING(11:80),fmt='(2F13.4)') FITLIM(1),FITLIM(2) C CALL sttput(STRING,STAT) do II = FITLIM(1),FITLIM(2) C STRING = ' II = ' C write (STRING(7:80),fmt='(i6)') II C CALL sttput(STRING,STAT) DX = II - XCONT(1) DY = DX * SLOPE GSUM = YCONT(1) + DY X0 = dble(II) C STRING = ' X0 = ' C write (STRING(7:80),fmt='(F13.4)') X0 C CALL sttput(STRING,STAT) do JJ = 1,NLINES AA(JJ) = dble(CENT(JJ)) AA(JJ+NLINES) = dble(SIGMA(JJ)) AA(JJ+2*NLINES) = dble(CMAX(JJ)) C ZZ = (X0 - CENT(JJ)) / SIGMA(JJ) C GFIT = CMAX(JJ) * exp(-(ZZ**2)/2) C GSUM = GSUM + GFIT end do if (MODE.ge.100) then C STRING = 'WRONG MODE' C CALL sttput(STRING,STAT) CALL lorentz(X0,Y0,AA,NLINES) else CALL gauss(X0,Y0,AA,NLINES) end if GSUM = GSUM + real(Y0) PIX = II-FITLIM(1)+1 C STRING = ' OUT( ) = ' C write (STRING(6:11),fmt='(i6)') PIX C write (STRING(16:80),fmt='(f13.4)') GSUM C CALL sttput(STRING,STAT) OUT(PIX) = GSUM end do return end SUBROUTINE SNLFIT(MODE,X,Y,W,N,NLINES,FUNC, & TT,T0,TSIG,D,WAV,VAR,S) C C X(N),Y(N) x,y-points I=1...N C W(N) = weights C Use MODE < 10 if W(I)=1, 10 < MODE < 20 if W(I) = 1/N (Poisson noise) C N = Number of points (less or equal 50000) C M = Number of parameters A(M) for fit function FUNC C (less or equal 40) C A(M) : start of SNLFIT: estimated initial parameters C end of SNLFIT: fitted parameters C A0(k): iteration limits for parameters A(k) C SIG(k) : Error estimation of parameters C FUNC : Name of Function C C H.-H. LOOSE US-GOETTINGEN Feb 1986 C changed 22 May 1987 C output lun=7 new 20 Jan 1990 C C This Version Mmax=40, Nmax=50000 C IMPLICIT DOUBLE PRECISION(A-H,O-Z), INTEGER(I-N) character*80 STRING !holder for output messages INTEGER STAT C DIMENSION X(1),Y(1),W(1) DIMENSION Y1(50000),Y2(40),A1(40),A2(40),D2(40),X2(2) DIMENSION A(40),A0(40),SIG(40),D(40),TT(40),T0(40),TSIG(40) DIMENSION G(40,41) EXTERNAL FUNC DATA ITMAX/100/ data IFLAG/0/ C C startup C STRING = ' SNLFIT Vers: 22-May-87' CALL STTPUT(STRING,STAT) C C initializations C if (MODE.ge.10) then IFLAG = 1 MODE = mod(MODE,10) end if if (MODE.eq.2) then M = 2*NLINES do II = 1,M A(II) = TT(II+NLINES) A0(II) = T0(II+NLINES) end do else if (MODE.eq.3) then M = 2*NLINES do II = 1,NLINES A(II) = TT(II) A0(II) = T0(II) end do do II = NLINES+1,M A(II) = TT(II+NLINES) A0(II) = T0(II+NLINES) end do else if (MODE.eq.4) then M = 2*NLINES + 1 do II = 1,NLINES+1 A(II) = TT(II) A0(II) = T0(II) end do do II = NLINES+2,M A(II) = TT(II+NLINES-1) A0(II) = T0(II+NLINES-1) end do else M = 3*NLINES do II = 1,M A(II) = TT(II) A0(II) = T0(II) end do end if ITER=0 DO 10 K=1,M 10 D2(K) = 2.D0*D(K) E=1.D-3 M1=M+1 C DN = N NFREE = N-M FREE = NFREE WAV = 0.D0 DO 11 I=1,N 11 WAV = WAV + W(I) WAV = WAV/DN C GOTO 2000 12 CONTINUE X2(1) = S C C - - main loop C 500 ITER = ITER + 1 STRING = ' Iteration ' WRITE(STRING(12:80),fmt='(i4)') ITER CALL STTPUT(STRING,STAT) C C more initializations C DO 30 K=1,M A1(K) = A(K) DO 30 J=K,M1 30 G(K,J) = 0.D0 C DO 40 I=1,N X0 = X(I) C - - derivatives DO 41 K=1,M A1(K) = A(K)*(1.D0+D(K)) CALL tmake(MODE,NLINES,A1,TT) CALL func(X0,Y3,TT,NLINES) A1(K) = A(K)*(1.D0-D(K)) CALL tmake(MODE,NLINES,A1,TT) CALL FUNC(X0,Y0,TT,NLINES) Y2(K) = (Y3-Y0)/(A(K)*D2(K)) A1(K) = A(K) 41 CONTINUE C - - matrix DO 42 K=1,M G(K,M1) = G(K,M1) + W(I)*(Y(I)-Y1(I))*Y2(K) DO 42 J=K,M 42 G(K,J) = G(K,J) + W(I)*Y2(K)*Y2(J) 40 CONTINUE DO 43 K=2,M DO 43 J=1,K-1 43 G(K,J) = G(J,K) C C CALL VMAT(G,M) C 1000 CONTINUE CALL SOLS(G,M,E,A2,SIG) C C update fit parameters C DO 50 K=1,M 50 A(K) = A1(K) + A2(K) I0 = 1 GOTO 2000 13 CONTINUE X2(2) = S IF( X2(2) .LT. X2(1) ) GOTO 14 E = E*1.D1 GOTO 1000 14 CONTINUE E = E/1.D1 X2(1) = X2(2) C L = 0 DO 60 K=1,M IF(DABS(A2(K)) .GT. DABS(A0(K))) L = 1 60 CONTINUE IF (ITER .GE. ITMAX) THEN STRING = ' *** Iteration limit ITMAX= reached ***' write (STRING(28:31),fmt='(i4)') ITMAX CALL STTPUT(STRING,STAT) L = 0 END IF IF(L .GT. 0) GOTO 500 C C - - end of main loop C STRING = ' Variance = Wmean=' write (STRING(12:24),fmt='(1PD12.4)') VAR write (STRING(35:80),fmt='(1PD12.4)') WAV CALL STTPUT(STRING,STAT) IF (IFLAG.eq.1) THEN FAC = 1.D0 ELSE FAC = VAR END IF DO 65 K=1,M 65 SIG(K) = FAC*DSQRT(SIG(K)) CALL tmake(MODE,NLINES,SIG,TSIG) CALL tmake(MODE,NLINES,A,TT) STRING = ' SNLFIT: END ITER=' write (STRING(20:80),fmt='(i4)') ITER CALL STTPUT(STRING,STAT) RETURN C C - - function values and chi**2 C 2000 S = 0 DO 70 I = 1,N X0 = X(I) CALL tmake(MODE,NLINES,A,TT) CALL FUNC(X0,Y0,TT,NLINES) Y1(I) = Y0 S = S + W(I)*(Y(I) - Y0)**2 70 CONTINUE IF(NFREE .GT. 0) S = S/FREE VAR = DSQRT(S/WAV) IF(I0 .EQ. 0) GOTO 12 GOTO 13 C END C ------------------------------------ subroutine vmat(g,m) implicit DOUBLE PRECISION(a-h,o-z), integer(i-n) dimension g(10,11) C m1 = m+1 do 100 iz=1,m write(6,300) (g(iz,is), is=1,m1) 300 format(20(f11.2)) 100 continue C return end subroutine sols(g,m,e,x,s) c c to be used with SNLFIT version 22-May-87 or later c c input: g(m,m+1) c m less or equal 40 c e Marquardt factor c output: x(m) solution c s(m) diagonal elements of inverse of g(m,m) c c H.-H. Loose US-Goettingen, 22-May-87 c implicit DOUBLE PRECISION(a-h,o-z), integer(i-n) dimension g(40,41), a(40,40), x(1), s(1) c m1 = m + 1 e1 = e + 1 c do 11 i=1,m do 10 k=1,m 10 a(i,k) = g(i,k) / dsqrt(g(i,i)*g(k,k)) 11 a(i,i) = e1 c CALL matinv(a,m,det) if (det .eq. 0.d0) stop 'SOLS' c do 20 j=1,m x(j) = 0.d0 s(j) = a(j,j) / g(j,j) do 20 k=1,m 20 x(j) = x(j) + g(k,m1)*a(j,k)/ dsqrt(g(j,j)*g(k,k)) c return end subroutine matinv (array, norder, det) c c purpose c invert a symmetric matrix and calculate its determinant c c description of parameters c array - input matrix which is replaced by its inverse c norder - degree of matrix (order of determinant) c det - determinant of input matrix c c subroutines and functions required c none c c comments c dimension statement valid for order up to 40 c c adapted from P.R. Bevington program B-2 c by H.-H. Loose checked and tested: 22-May-1987 c implicit DOUBLE PRECISION (a-h,o-z), integer (i-n) dimension array(40,40), ik(40), jk(40) 10 det = 1.d0 11 do 100 k=1,norder c c find largest element array(i,j) in rest of matrix c amax = 0.d0 21 do 30 i=k,norder do 30 j=k,norder 23 if(dabs(amax)-dabs(array(i,j))) 24,24,30 24 amax = array(i,j) ik(k) = i jk(k) = j 30 continue c c interchange rows and columns to put amax in array(k,k) c 31 if (amax) 41,32,41 32 det = 0.d0 goto 140 41 i = ik(k) if (i-k) 21,51,43 43 do 50 j=1,norder save = array(k,j) array(k,j)=array(i,j) 50 array(i,j) = -save 51 j = jk(k) if (j-k) 21,61,53 53 do 60 i=1,norder save = array(i,k) array(i,k) = array(i,j) 60 array(i,j) = -save c c accumulate elements of inverse matrix c 61 do 70 i=1,norder if (i-k) 63,70,63 63 array(i,k) = -array(i,k)/amax 70 continue 71 do 80 i=1,norder do 80 j=1,norder if (i-k) 74,80,74 74 if (j-k) 75,80,75 75 array(i,j) = array(i,j) + array(i,k)*array(k,j) 80 continue 81 do 90 j=1,norder if (j-k) 83,90,83 83 array(k,j) = array(k,j)/amax 90 continue array(k,k) = 1.d0/amax 100 det = det*amax c c restore ordering of matrix c 101 do 130 l=1,norder k = norder - l + 1 j = ik(k) if (j-k) 111,111,105 105 do 110 i=1,norder save = array(i,k) array(i,k) = -array(i,j) 110 array(i,j) = save 111 i = jk(k) if (i-k) 130,130,113 113 do 120 j=1,norder save = array(k,j) array(k,j) = -array(i,j) 120 array(i,j) = save 130 continue 140 return end subroutine GAUSS(X0,Y0,AA,NLINES) implicit none DOUBLE PRECISION X0,Y0 !data values DOUBLE PRECISION AA(40) !fit parameters INTEGER NLINES !number of gaussian profiles being fit DOUBLE PRECISION CENT !center of line DOUBLE PRECISION SIGMA !gaussian sigma of line DOUBLE PRECISION CMAX !line maximum INTEGER KK !dummy counter DOUBLE PRECISION ZZ,GFIT !dummy value holder Y0 = 0. do KK = 1,NLINES CENT = AA(KK) SIGMA = AA(KK+NLINES) CMAX = AA(KK+2*NLINES) if (SIGMA.ne.0) then ZZ = (X0 - CENT) / SIGMA C Cut-off at GFIT < 1e-6 * CMAX C if (ZZ.gt.5.3) then C GFIT = 0 C else GFIT = CMAX * exp(-(ZZ*ZZ)/2) C end if Y0 = Y0 + GFIT end if end do return end subroutine LORENTZ(X0,Y0,AA,NLINES) implicit none DOUBLE PRECISION AA(40) !fit parameters DOUBLE PRECISION BETA2 !(FWHM/2)**2 DOUBLE PRECISION CENT !center of line DOUBLE PRECISION CMAX !line maximum DOUBLE PRECISION FWHM !gaussian sigma of line INTEGER NLINES !number of gaussian profiles being fit DOUBLE PRECISION X0,Y0 !data values INTEGER II !dummy counter DOUBLE PRECISION ZZ,GFIT !dummy value holders Y0 = 0. do II = 1,NLINES CENT = AA(II) FWHM = AA(II+NLINES) CMAX = AA(II+2*NLINES) BETA2 = (FWHM / 2)**2 ZZ = BETA2 / ((X0 - CENT)**2 + BETA2) GFIT = CMAX * ZZ Y0 = Y0 + GFIT end do return end subroutine TMAKE(MODE,NLINES,AA,TT) C reconstruct TT matrix of params from AA (possibly partial) matrix implicit none DOUBLE PRECISION AA(40) !fit parameters INTEGER MODE !fit mode INTEGER NLINES !number of gaussian profiles being fit DOUBLE PRECISION TT(40) !function parameters INTEGER MM !number of parameters in AA INTEGER II !dummy counter if (MODE.eq.2) then MM = 2*NLINES do II = 1,MM TT(II+NLINES) = AA(II) end do else if (MODE.eq.3) then MM = 2*NLINES do II = 1,NLINES TT(II) = AA(II) end do do II = NLINES+1,MM TT(II+NLINES) = AA(II) end do else if (MODE.eq.4) then MM = 2*NLINES + 1 do II = 1,NLINES+1 TT(II) = AA(II) end do do II = NLINES+2,2*NLINES TT(II) = AA(NLINES+1) end do do II = NLINES+2,MM TT(II+NLINES-1) = AA(II) end do else MM = 3*NLINES do II = 1,MM TT(II) = AA(II) end do end if return end