C @(#)fitlyman.for 10.4 (ESO-IPG) 2/9/96 12:10:18 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 fitlyman implicit none include 'MID_REL_INCL:fit_var.inc' ! copies of keywords coming from line character*80 SPESYS,SESSYS,OUTSYS,HISSYS integer GRASYS character*30 choice character*70 sjunk,ErrMes character*60 filnam integer ifla,ist,i,ijk,icur,IDold,iok real xdum,ydum,zdum,vdum,c parameter (c=2.997e5) logical la C include 'MID_REL_INCL:fit_mid.inc' INTEGER MADRID(1) INCLUDE 'MID_INCLUDE:ST_DEF.INC' COMMON/VMR/MADRID INCLUDE 'MID_INCLUDE:ST_DAT.INC' C!!!!!!! C MIDAS SET-UP C!!!!!!! call STSPRO('fitlyman') call STKRDC('SPECTAB',1,1,80,i,SPESYS,ijk,ijk,iok) if (iok.ne.0) then call ErrMsg('Error in reading keywords: check context..') call STSEPI() stop end if call STKRDC('OUTTAB',1,1,80,i,OUTSYS,ijk,ijk,iok) if (iok.ne.0) then call ErrMsg('Error in reading keywords: check context..') call STSEPI() stop end if call STKRDC('HISTAB',1,1,80,i,HISSYS,ijk,ijk,iok) if (iok.ne.0) then call ErrMsg('Error in reading keywords: check context..') call STSEPI() stop end if call STKRDC('SESSNAM',1,1,80,i,SESSYS,ijk,ijk,iok) if (iok.ne.0) then call ErrMsg('Error in reading keywords: check context..') call STSEPI() stop end if call STKRDI('GRAWIN',1,1,i,GRASYS,ijk,ijk,iok) if (iok.ne.0) then call ErrMsg('Error in reading keywords: check context..') call STSEPI() stop end if C!!!!!!! C Videata iniziale C!!!!!!! choice = ' ' do i=1,4 call sttdis(choice,0,ist) end do sjunk = ' - FitLyman- Lyman clouds fitting procedure' do i=1,2 call sttdis(' ',0,ist) end do call DisMsg('Loading & initializing data...') call DisMsg('Wait please..') C!!!!!!! C Inizializzazione programma C!!!!!!! if (SESSYS.ne.'NULL') then i=index(SESSYS,' ') sjunk=SESSYS(1:i-1)//PARTBL else sjunk=PARTBL end if call rsetup(sjunk,ist) ! read set-up if (ist.ne.0) then ! on error, default values call WrnMsg('Error in reading set-up: using defaults') ! copy syskeys into variables if (SPESYS.ne.'NULL') then FILSPE=SPESYS else call ErrMsg('Missing input spectrum file') call STSEPI() stop end if if (OUTSYS.ne.'NULL') then FILOUT=OUTSYS else i=index(FILSPE,' ') FILOUT=FILSPE(1:i-1)//'FIT' end if if (HISSYS.ne.'NULL') then FILLOG=HISSYS else FILLOG=FILSPE end if TURBLN=0 i_grap = .false. if (GRASYS.ne.0) i_grap = .true. call F_TBLR('SCRATCH',1,ist) !initializes lypar (in memory!) else call f_tblr(sjunk,1,ist) !read init. param. if (SPESYS.ne.'NULL') FILSPE=SPESYS ! overrides values read if (OUTSYS.ne.'NULL') FILOUT=OUTSYS if (HISSYS.ne.'NULL') FILLOG=HISSYS if (GRASYS.ne.-1) then if (GRASYS .eq.1) i_grap = .true. if (GRASYS .eq.0) i_grap = .false. end if end if call f_tblw(PARTBL,1,ist) c call AskStp('PROGRAM',ist) call ssetup(PARTBL,ist) if (ist.ne.0) then call ErrMsg('Error in writing set-up: check disk space') call STSEPI() stop end if call AtmRD(ist) !read atompar.dat NPUNTI=NMXSPE call reaspe(FILSPE,NPUNTI) !read spectrum if (NPUNTI.lt.1) then call STSEPI() stop endif i=index(FILLOG,' ') filnam=fillog(1:i-1)//PARTBL call GETIDN(filnam,IDnum,ist) if (ist.ne.0) then call WrnMsg('Log files not found: starting from scratch') IDnum=0 end if ! read fit limits if (SESSYS.ne.'NULL') then i=index(SESSYS,' ') sjunk=SESSYS(1:i-1)//LIMTBL else sjunk=LIMTBL end if call rintvl(sjunk,1,ist) call sintvl(sjunk,1,ist) ! read minuit instructions if (SESSYS.ne.'NULL') then i=index(SESSYS,' ') sjunk=SESSYS(1:i-1)//MINTBL else sjunk=MINTBL end if call rminui(sjunk,1,ist) call sminui(sjunk,1,ist) call ReaRes call GraMai(ist) ! plot call DisMsg(' ..done') call sttdis(' ',0,ist) C!!!!!!! C Menu pricipale C!!!!!!! 1001 continue choice = ' ' call MMenu(choice) C#### C SAVESESSION C#### if (choice.eq.'SAVESESSION') then ! saves a session ifla=-2 call AskC('Session name',sjunk,ifla) i=index(sjunk,' ') filnam=sjunk(1:i-1)//PARTBL call f_tblw(filnam,1,ifla) call ssetup(filnam,ifla) filnam=sjunk(1:i-1)//LIMTBL call sintvl(filnam,1,ifla) filnam=sjunk(1:i-1)//MINTBL call sminui(filnam,1,ifla) end if C#### C RECOVER C#### if (choice.eq.'RECOVER') then ! load a session ifla=-2 call AskC('Session name',sjunk,ifla) i=index(sjunk,' ') filnam=sjunk(1:i-1)//PARTBL call f_tblr(filnam,1,ifla) call rsetup(filnam,ifla) filnam=sjunk(1:i-1)//LIMTBL call rintvl(filnam,1,ifla) filnam=sjunk(1:i-1)//MINTBL call rminui(filnam,1,ifla) call GraMai(ifla) call ShPar1('NORMAL') call f_tblw(PARTBL,1,ifla) call ssetup(PARTBL,ifla) ! saves on std tbl call sintvl(LIMTBL,1,ifla) call sminui(MINTBL,1,ifla) end if C#### C SETUP C#### if (choice.eq.'SET-UP') then ! ASK FOR SETUP VALUES call AskStp(' ',ifla) if (ifla.eq.0) then call ssetup(PARTBL,ifla) !saves new setup else call rsetup(PARTBL,ifla) ! load old setup end if call GraMai(ifla) end if C#### C END C#### if (choice.eq.'END') then if (I_GRAP) then call PTCLOS () end if c close(IHFILE_LOG) call STSEPI() stop end if C#### C NEW C#### if (choice.eq.'NEWLINE') then ! load a scratch session call f_tblr('SCRATCH',1,ifla) call rintvl('SCRATCH',1,ifla) call rminui('SCRATCH',1,ifla) call f_tblw(PARTBL,1,ifla) ! saves on std tbl call ssetup(PARTBL,ist) call sintvl(filnam,1,ifla) call sminui(filnam,1,ifla) end if C#### C ITERATE C#### if (choice.eq.'ITERATE') then ! load from results call ShoRes ! Carica l'ultimo fit do i=1,NROWS LamIni(i)=LamCen(i) NIni(i)=NHfit(i) BIni(i)=BDopp(i) BTuIni(i)=BTur(i) end do call ShPar1('NORMAL') call f_tblw(PARTBL,1,ifla) ! saves on std tbl call ssetup(PARTBL,ist) call sintvl(LIMTBL,1,ifla) call sminui(MINTBL,1,ifla) end if C#### C HISTORY C#### if (choice.eq.'HISTORY') then ! load from logfiles call AskI('Input ID number to be loaded:',IDold,iok) if (iok.ne.-1) then i=index(FILLOG,' ') filnam=FILLOG(1:i-1)//PARTBL !reads from logfile call f_tblr(filnam,IDold,ist) if (ist.ne.0) then call ErrMsg('Error in reading parameter table') call f_tblr(PARTBL,1,ist) ! reads standard tables call rminui(MINTBL,1,ist) call rintvl(LIMTBL,1,ist) goto 1001 end if filnam=FILLOG(1:i-1)//MINTBL !reads from logfile call rminui(filnam,IDold,ist) if (ist.ne.0) then call ErrMsg('Error in reading minuit command table') call f_tblr(PARTBL,1,ist) ! reads standard tables call rminui(MINTBL,1,ist) call rintvl(LIMTBL,1,ist) goto 1001 end if filnam=FILLOG(1:i-1)//LIMTBL !reads from logfile call rintvl(filnam,IDold,ist) if (ist.ne.0) then call ErrMsg('Error in reading interval table') call f_tblr(PARTBL,1,ist) ! reads standard tables call rminui(MINTBL,1,ist) call rintvl(LIMTBL,1,ist) goto 1001 end if call ShPar1('NORMAL') !Shows new parameters call f_tblw(PARTBL,1,ist) ! save on standard tables call ssetup(PARTBL,ist) call sminui(MINTBL,1,ist) call sintvl(LIMTBL,1,ist) end if end if C#### C STANDARD C#### if (choice.eq.'STANDARD') then 501 call AskPar(ifla) if (ifla.ne.0) then !quit call f_tblr(PARTBL,1,ifla) if (ifla.ne.0) call ErrMsg('Error in reading lypar.tbl') goto 1001 end if call ChkPar(ErrMes,ifla) if (ifla.ne.0) then call ErrMsg(ErrMes) goto 1001 end if call AddAtP(ifla) if (ifla.ne.0) then call ErrMsg('Atomic parameter not found') goto 1001 end if call AdCoef() call AddInp() call AdPLim() call f_tblw(PARTBL,1,ifla) if (ifla.ne.0) then call ErrMsg('Error in writing lypar.tbl') goto 1001 end if call ssetup(PARTBL,ist) 511 call AskFIn(ifla) if (ifla.ne.0) then !quit call rintvl(LIMTBL,1,ifla) goto 501 else call sintvl(LIMTBL,1,ifla) end if call AskMin(ifla) if (ifla.ne.0) then !quit call rminui(MINTBL,1,ifla) goto 1001 !goes back to main menu else call sminui(MINTBL,1,ifla) end if call sttdis(' ',0,ifla) call DisMsg('Preparing minimization data...') call FitWin(i,ifla) ! create data for FCN if (ifla.lt.0) then call ErrMsg('Error in initial condition for fit') goto 1001 ! torna al menu' principale end if call FCNHd() call MinuHd() call DisMsg('...Done') write(sjunk,'(a25,i6)')'Number of pixel selected:',i call DisMsg(sjunk) C Esegue il fit call DisMsg('Executing fitting procedures...') call Minmze() call DisMsg('...Done') C Presenta e salva i risultati call ShoRes call GraMai(ifla) C SAVE? call sttdis(' ',0,ist) call AskYN('Save results in work & log files? [Y/N]', 1 la,ist) if (la) then IDnum=IDNum+1 call SavRes(FILOUT,ist) i=index(FILLOG,' ') !saves on history file filnam=FILLOG(1:i-1)//PARTBL call f_tblw(filnam,IDNum,ist) filnam=FILLOG(1:i-1)//MINTBL call sminui(filnam,IDNum,ist) filnam=FILLOG(1:i-1)//LIMTBL call sintvl(filnam,IDNum,ist) write(sjunk,'(a17,i6)')'Saved with ID = ',IDnum call DisMsg(sjunk) end if end if C#### C EDITPARAM C#### if (choice.eq.'EDITPARAM') then 601 call EdtPar(ifla) if (ifla.ne.0) then !quit call f_tblr(PARTBL,1,ifla) if (ifla.ne.0) then call ErrMsg('Error in reading lypar') goto 1001 end if else call sttdis(' ** Checking consistency of lypar.tbl...', 1 0,ist) call ChkPar(ErrMes,ifla) if (ifla.ne.0) then call ErrMsg(ErrMes) goto 601 ! edit again.. end if call sttdis(' ...done' 1 ,0,ist) call AddAtP(ifla) if (ifla.ne.0) then call ErrMsg('Atomic parameter not found') goto 601 ! edit again.. end if call AdCoef() call AddInP() call AdPLim() call f_tblw(PARTBL,1,ifla) if (ifla.ne.0) then call ErrMsg('Error in writing lypar') goto 1001 end if call ssetup(PARTBL,ist) end if end if C#### C EDITLIM C#### if (choice.eq.'EDITLIM') then call AskFIn(ifla) if (ifla.ne.0) then call rintvl(LIMTBL,1,ifla) else call sintvl(LIMTBL,1,ifla) end if end if C#### C EDITMINUIT C#### if (choice.eq.'EDITMINUIT') then call AskMin(ifla) if (ifla.ne.0) then call rminui(MINTBL,1,ifla) else call sminui(MINTBL,1,ifla) end if end if C#### C DIRECTMINI C#### if (choice.eq.'DIRECTMINI') then call DisMsg('Preparing minimization data...') call f_tblr(PARTBL,1,ifla) if (ifla.ne.0) then call ErrMsg('Error in reading lypar') goto 1001 end if call rintvl(LIMTBL,1,ifla) if (ifla.ne.0) then call ErrMsg('Error in reading lypar') goto 1001 end if call rminui(MINTBL,1,ifla) if (ifla.ne.0) then call ErrMsg('Error in reading lypar') goto 1001 end if call FitWin(i,ifla) ! create data for FCN if (ifla.lt.0) then call ErrMsg('Error in initial condition for fit') goto 1001 ! torna al menu' principale end if call FCNHd() call MinuHd() call DisMsg('...Done') write(sjunk,'(a25,i6)')'Number of pixel selected:',i call DisMsg(sjunk) C Esegue il fit call DisMsg('Executing fitting procedures...') call Minmze() call DisMsg('...Done') C Presenta e salva i risultati call ShoRes call GraMai(ifla) C SAVE? call sttdis(' ',0,ist) call AskYN('Save results in work & log files? [Y/N]', 1 la,ist) if (la) then IDnum=IDNum+1 call SavRes(FILOUT,ist) i=index(FILLOG,' ') !saves on history file filnam=FILLOG(1:i-1)//PARTBL call f_tblw(filnam,IDNum,ist) filnam=FILLOG(1:i-1)//MINTBL call sminui(filnam,IDNum,ist) filnam=FILLOG(1:i-1)//LIMTBL call sintvl(filnam,IDNum,ist) write(sjunk,'(a17,i6)')'Saved with ID = ',IDnum call DisMsg(sjunk) end if end if C#### C DEFINEWINDOW C#### if (choice.eq.'DEFINEWINDOW') then call DefGrW(ifla) if (ifla.eq.0) then call ssetup(PARTBL,ifla) call GraMai(ifla) end if endif C#### C CURSOR C#### if (choice.eq.'CURSOR') then call GetXCr(xdum,ydum,zdum,vdum,icur) end if C#### C GOBACK C#### if (choice.eq.'GOBACK') then if (I_Z) then REDSH=(-2.*VELRAN*4./5.-REDSH*VELRAN*4./5.+2.*c*REDSH) 1 /(2.*c + VELRAN*4./5.) else do ijk = 1,NREGIO REGMIN(ijk) = REGMIN(ijk) - WINSTE(ijk)*4./5. REGMAX(ijk) = REGMAX(ijk) - WINSTE(ijk)*4./5. end do end if call GraMai(ifla) call ssetup(PARTBL,ifla) !saves new setup endif C#### C GOFORW C#### if (choice.eq.'GOFORW') then if (I_Z) then REDSH=(2.*VELRAN*4./5.+REDSH*VELRAN*4./5.+2.*c*REDSH) 1 /(2.*c - VELRAN*4./5.) else do ijk = 1,NREGIO REGMIN(ijk) = REGMIN(ijk) + WINSTE(ijk)*4./5. REGMAX(ijk) = REGMAX(ijk) + WINSTE(ijk)*4./5. end do end if call GraMai(ifla) call ssetup(PARTBL,ifla) !saves new setup end if goto 1001 end