SPECTRUM - bug reports and modification history 2001 -- C.S.Jeffery --------------------------------------------------------------------------- Summary as of 2003.06.01 Major changes: - SPECTRUM v3 introduced - DP_LINES can preselect lines from linelist - SP_SYNTH can treat up to 500000 frequencies. - depth-dependent velocities can be considered. - LSYNT is now much more efficient. Major bugs fixed: - LHI: microturbulence > 0 is now working - Ca H+K are now correctly treated - ** anything heavier than iron was WRONG before 22 Nov 2001. fixed. ** Persistent problems: - CONFL: overflow in C/L ratio - ATSCAN: requires check for existing ATSCAN Matters to attend to: - HEI, HEII: do not include microturbulence - H does not include microturbulence in synthesis - LMETAL, LSYNT: (b) van der Waal's broadening incomplete - HEII: Need to include all of the Schoening Butler tables. Long-term plans: - ATSCAN: Preselect lines to include in synthesis on basis of predicted line strength (in order to save time computing very weak lines) - ATLAS_RD, TLUSTY_RD: Include options to read "foreign" models directly - ATSCAN: should set default abundances to those of model atmosphere. - Port to Linux has still to be completed ------------------------------------------------------------------------ New journal commenced 2001.Nov.15. Please place new entries at head of file. ------------------------------------------------------------------------ 2003 June 27 LMETAL/LHI/LHEI/LHEII Located references to "U" and replaced with TAP_PARFN. Tested LMETAL, but not the rest. u.f in library dp should now be redundant. SPC_STARTUP Aaagh.... Where have the HeI lines gone ???? Looks like the Beachamp tables were not being read with nhedat = 2. HeI lines OK for BCS tables. -- FIXED - must read in both sets of data at startup - so that nhedat can be changed at will. 2003 Jun 26 U This subroutine should have been replaced by TAP_PARFN years ago. However there are still several calls in various subroutines. I have replaced U with a dummy front end to TAP_PARFN, which should fix the problems. Note: there are still calls to U to replace. -- allegedly done -- 02/07/2003 However Fe abundances now different to old version, and old version was preferred. Problem was ZION - needs to be 0 for neutrals. -- fixed and sent for testing -- 02/07/2003 EWABUND Causing problems in some casese with poor convergence in SPECTRUM v3. May be related to problems in INTEG_FEAUT. Divergence trapped and flagged. Code and logic now need tidying up. INTEG_FEAUT Single line mode for some lines produces broad "emission" wings. I assume that this is due to poor treatment of the continuum opacity in this mode. It needs to be corrected. -- done -- 02/07/2003 -- WHOLE The "calc_abund" option needs to use a different version of WHOLE than that used by other modes. Currently done by renaming the new (internal) WHOLE as SPC_WHOLE. The original WHOLE will then be called as before. 2003 Jun 3 DP_HGRID replaces HGRID and HSEARCH. The latter failed to include all the high-order series members in a synthesis because of a terrible algorithm. 2003 May 21 - 31 DP_LINES Implemented line selection procedure. Uses a parameter line_threshold to eliminate weak lines SP_SCAN Introduced HEAPSORT to sort lines in long lists, because existing sort very inefficient for NELSY > 50. THIS IS NOT YET OPERATIONAL. SPECTRUM various other items implemented to deal with depth-dependent velocity fields. 2003 May 16 - 20 : SPECTRUM 3 TESTS : /electra/models/test/spectrum3/Test.sh 1) Mode 0: Balmer Line profiles Teff=20000, log g = 4.0, nH=12.00: OK 2) Mode 0: Balmer Line abundances Teff=20000, log g = 4.0, EW=...: OK 3) Mode 10: Spectrum synthesis standard abundances Teff=20000, log g = 4.0, WRANGE={4300., 4500.), vt=0.0: The two versions produce slightly different wavelength grids - by default, for whatever reason. The difference in integrated equivalent width over 200 Angstroms was less than 0.07 per cent. 4) Mode 11: Spectrum synthesis modified abundances Teff=20000, log g = 4.0, WRANGE={4300., 4500.), vt=5.0: H = 11.7, He = 11.3, C = 9.3 Perfect agreement. 5) Mode 11: As (4) with specific intensities Teff=20000, log g = 4.0, WRANGE={4300., 4500.), vt=5.0: H = 11.7, He = 11.3, C = 9.3 Perfect agreement. Note that for mu=0.005, the intensities do not look very good for EITHER new or old calculations. This is probably because the flux integral calculation breaks down badly when the optical depth changes rapidly over zones -- ie tau_line is undersampled. ------------------------------------------------------------------------ 2003 May 15: SPECTRUM 3 I have embarked on a major rewrite of Spectrum. The following principal features are introduced. i. Some degree of backward compatability is maintained, meaning that old input scripts with prescribed numeric only inputs should still work. These are obsolescent, and cannot be guaranteed for a future release of SPECTRUM. ii. The code will be generally driven by a command structure, thus flags and switches are set specifically in the input sequence, and a dictionary of commands and keywords will be provided. iii. In order to provide greater flexibility, the subroutine WHOLE has been moved to become an internal subroutine of the SPECTRUM. This allows more control without the introduction of global variables. iv. In future, there will be more rationalisation of global variables, through COMMON, and a phased transition to use more features of fortran90. v. These changes have been made in order, partly, to be able to introduce depth-dependent quantities in the the models in a logical way, rather than by making ad hoc or obscure tricks in the input code. vi. The old "prescribed input" sections of SPECTRUM are now an internal subroutine of SPECTRUM, called SPC__V1. The selection of which version is used is still made by the initial input line. i1 i2 r3 i4 If i1 > 0, old input is assumed If i1 = 0, new style input is assumed. WHOLE subroutine has been moved to become an internal subroutine of SPECTRUM. DP_VDATA a new subroutine provided to map a prescribed velocity structure (either radial or mirocturbulent) given as a function of optical depth onto the model atmopshere optical depth structure. Global arrays QUB__VRAD and QUB__VTRB contain the result. OUTDIP Modified to allow overwriting of previous dipso.xxx files. DP_ATDATA split up so that DP_ATDATA reads in general atomic data, ionisation energies, cosmic abundances and so on. DP_LINES new subroutine to read in the linelist. Now, we can introduce a facility to read in the model first, then read in the linelist and preselect lines as they are read in on the basis of strength. Not yet implemented but to do soon. Takes an arguemnt for the filename. MODEL_ATM will now take a model filename as argument. If fuinemae is an empty string, then subroutine is backward compatabile; file must be opened before first call, and is not closed at end of read. If called with an argument, expects only one model per file. SP_SYNTH modified to take above structural changes into account. 2003 Mar 21 DP_RD_IONFRN If a file IONFRN exists in the current directory, MODEL_ATM will be directed to read the ionisation fractions and partition functions from this file. The purpose is to allow non-LTE model atmospheres to be used in the line formation calculation. i.e. the non-LTE ionisation fractions (and self-consistent partition functions) are used, although within an ion, the level populations will still follow the Saha-Boltzmann values. See the subroutine for the format of the file. 2003 Mar 15: lines_woolf.d Vincent Woolf has supplied a new compilation of atomic lines which includes all lines in the blue-optical likely to contribute to the absorption spectrum of extreme hlium stars. He has described the procedure he used to construct this list below. This file (lines_woolf.d) is packaged with the distribution of SPECTRUM. -- To create it I compiled all the lines in the Kurucz list between 3800 and 5200 A. I then ran them through Spectrum for twelve Teff, logg settings to cover the range likely for EHe's. Model atmospheres were from /electra/models/009901 (higher metals). parameters for tests: T/1000 logg 10 0.75 10 3.0 10 4.0 12 4.0 15 3.5 20 1.5 20 3.5 20 4.0 24 3.5 26 4.0 30 4.0 30 2.5 Lines which produced equivalent widths greater than 1.0 mA in any test run were included in the line list. It's possible that groups of weaker than 1.0 mA lines could combine to make important features in a spectrum, but it'll be a rare event and there had to be some cutoff to reduce the number of lines. I made some effort to combine lines with the same wavelength, element, and excitation into single lines, adding the gf's. Might have missed a few, but this just means the number of lines in the list is larger, not that the spectrum output will be different. I then added in lines previously used in the Spectrum analyses (A2.d, Al2.d, Al3.d, C1.d, etc.) checking for lines duplicated between the csj and Kurucz lists. For some lines (fewer than 50) it was difficult to tell if there was a duplication. I made my best guess for these, but might have deleted a few lines that weren't duplicates or included a few that were. There are few enough like this that it probably won't make any noticible difference in the spectrum fits. But for analyses of individual lines it's a good idea to go back to the original, uncombined lists for the atomic data. I didn't preserve the reference information from the csj line lists. The references for the Kurucz data is on his web page, but not in his line list files. The "0.00 AAA 11/ /" bit is included to keep any formatted input routine happy. -- 2003 Mar 08: TAP_PARFN The old version of the partitiion function routine did not contain all of the Baschek et al. tables. This cause a problem when computing UV spectra for MnIII and CrIII. The subroutine has been updated and now contains ALL the Baschek et al. data. In addition, SPECTRUM was switched over to use the tap95/ (f95) library, although it still uses the old f77 TAP include files. 2003 Feb 27: QUB/OPHI I found a discrepancy between the Balmer edges given by the OP and those computed for standard model atmospherers. The Rydberg constant used in OPHI seems to be wrong a by a few per cent. It seems that OPK2/OPKHI did not have the same problem. Since OPHI only used teh Rydberg constant iin computing cutoff values for summing the opacity from each continuum I think it can be safely changed. Doing so corrected the problem. 2003 Feb 27: INTEG_FEAUT Continuum flux wavelength grid. a) This is now spaced in resolution increments according to wavelength so that UV spectra are better sampled than optical (per Angstrom). b) I have included a check that searches for opacity edges in a huge list of ionization ergies obtained from the Opacity Project. Any edges are mapped by placing a wavelength point either side of it. 2003 Feb 26: lte-codes Chris Winter has converted SPECTRUM + STERNE + SFIT to use the automake facility, which is much better for constructing a distribution version of the codes. This is now in the test stage. 2003 Feb 20: INTEG_FEAUT Changed from cubic Spline interpolation to linear interpolation for computing the local continuum fluxes. May be temporary until we can sort out what happens at opacity edges. 2003 Feb 20: MTH_D_TABINT New subroutine in mth library - double precision version of TABINT which performed newton poynomial interpolation of order (n) up to 10. 2003 Feb 19: Test Result: Model T25G40.d Linelist Default. Synthesis from 3300 -- 4000 => Three significant problems ---- 1. Continuum sampled every 10 A + spline interpolation led to bad oscillations at opacity edges. a) replace with linear interpolation as a temporary measure - DONE b) introduce extra wavelength points near opacity edges NB: INTEG_FEAUT computes continuum fluxes every 10A, but still computes the continuous opacity at every wavelength point. This could be better done by interpolation if we can store the continuous opacities first. 2. The Balmer series limit is not correctly reproduced. The high-order H lines are simply not in the code anywhere. Must fix. 3. A HUGE line at 3641.2 - error in the linelist (log gf = 7 !!): corrected . 2002 Dec 6: DP_ATDATA, ATABUND DP_ATDATA was not picking up element ids in TAP_ENAME because the intel fortran compiler was not setting a valid value for STATUS. This proagated through. Now fixed. Also made dialogue in ATABUND more economical. 2002 Jun 12: 2002 Jun 12: VCS Another dimension argarray(1) bug removed 2002 Apr 25: Port to Linux using the Intel compiler. Runs and produces realistic looking spectra. Head to head comparison produced a perfect match for H-rich models. VCS2, HEI_EWID, DP_ATDATA, TRANSFER, OUTSP ifc -cB threw up an array bounds violation . traced generally to declarations of array arguments in the obsolete form dimension argarray(1) 2002 Jan 14: TAP: A fortran 90 version of the TAP library has been released and tested with Spectrum. Works well. Note, however, that the "INCLUDE TAP..." statements in SPECTRUM require the F77 versions in $TOP/include COMMON_OPAC & COMMON_ALL: I have modified the include files QUB_CONTROL and created a new include file QUB_IONS so that COMMON_OPAC now points _only_ to other include files. COMMON_ALL points to the SAME inlcude files so that you whould not get problems from just editing one common in one inlcude file and forgetting to update the other. To make changes in a common block and then rebuild the program: >> cd $TOP/spectrum/common >> < do the edits > >> make install # this ensures the links are in ../include/ >> cd $TOP/pack/dp; make >> cd $TOP/pack/qub; make >> cd $TOP/spectrum; make 2001 Nov 28: Suggestions for changes: COMMON: needs to avoid duplicated common blocks in different INCLUDE files U & QAS: These subroutines should be removed and all references replaced either by references to U_ION, F_ION or TAP_PARFN. Subroutines affected are: LHEI, LHEI_B, LHEII, LHI. 2001 Nov 28: a) After major changes, the chi_lupi test stopped working. Finally traced it to a failure in the QUB opacity routines because I had changed some common blocks :( . b) There is also a problem in H line profiles at low Teff. The wings have a series of steps. Meanwhile I made incremental improvements to: DP/ DP_ATDATA (used to be ATDATA) DP_ATHEI (used to be ATHEI) DP_ATHEII (used to be ATHEII) LSYNT DP_ION DP_IONFRN OPAC (passes depth pointer to OPAQUB) VCS (extrapolation prevented: but did not solve "b". VCS2 (ditto). QUB/ OPAQUB (now uses pretabulated partition functions for Si and Mg) TAP/ TAP_IONPT COMMON/ QUB_DIMS (new file to contain dimensions for QUB and DP libraries) COMMON_ALL (includes QUB_DIMS) COMMON_OPAC (includes QUB_DIMS and DP__IONS ) 2001 Nov 23: DP_VDWC: New subroutine, called by MODEL_ATM Simply calculates depth-dependent part of van der Waal's broadening coefficient. Default values are "allegedly" computed in ATDATA. The results should be appied in LSYNT and LMETAL. They are currently not, because I do not trust/understand any of the formulae I am using. The code is simply put in place to ensure that adequate holders are present in common blocks, etc.... >>> Van der Waal's broadening is NOT YET IMPLEMENTED <<< LSYNT This now uses pretabulated functions; which is a great improvement. The old code used function "U", which returned a dummy value for Z>26! and was WRONG. The new version uses TAP_PARFN. This is VERY VERY important. TAP_PARFN There was an error in the Strontium data: ! Strontium. DATA (NS(J,38),J=1,4) / 2, 1, 1, 1 / === DP_ION: New subroutine, called by DP_IONFRN Simulates OPK_ION but is valid for more ions (Z>56!) and uses ionization potentials read in from TAP by ATDATA ATDATA Ensure a valid i.p. is available for all states up to 6. DO J = 1,6 RION(J,I) = TAP_IONPT ( I, J-1 ) !. if i.p. is not known (returned as zero) add 10 to previous i.p. IF (RION(J,I).LT.-999.D+00) RION(J,I) = RION(J-1,I)+10. ENDDO Default values for Gamma_s and Gamma_v are also provided. The former is probably OK (classical values from the impact approximation). The second needs more work. 2001 Nov 22: Following changes made: have to be tested tomorrow because I'm running a big grid today. Changes will be worthwhile because they will save _days_ of computing time for big grids. DP_IONFRN: New subroutine, called by MODEL_ATM Initializes array of ionization fractions and partition functions for all elements (currently up to Barium) to be used in LSYNT. Derived from OPK_IONFRN MODEL_ATM Calls DP_IONFRN COMMON_ALL 1) Includes ionization fraction arrays for DP_IONFRN 2) KA = 5 (formerly 100) in prder to save memory. KI = 50000 (formerly 10000) in order to include more lines 3) COMMON/IONS/ COSMIC(KI),ATW(KI),RION(6,KI),PARFUN(5,KI),IION(KI) ==== to allow for higher ionization species LSYNT New code to use F_ION rather than compute ion fraction for every line... Can later use this in ATDATA or ATSCAN to eliminate unwanted lines from linelist at an earlier stage and save even more CPU time. 2001 Nov 21: LSYNT Further improvements to optimize accuracy and efficiency. This involves ensuring that the line is integrated over a sufficient wavelenth interval, but not over too wide an interval (DOPCUT) and not if the line is insufficiently strong (GNFCUT). ! Wavelength cutoff in Doppler units for calculating metal lines ! Cutoff is DOPCUT*AV0*DOPP = DOPCUT*GAMMA(FWHM) Angstroms ! from line centre REAL DOPCUT PARAMETER (DOPCUT = 1.E+5) ! Cut off for n.gf to be used in line opacity calculation. ! the cutoff used is 10 times lower than would give a ! 0.4% error in the profiles of very weak lines. REAL GFNCUT PARAMETER (GFNCUT = 1.E-6) !.............................. ! ,only treat line if strong enough: IF (FRAC.GT.GFNCUT) THEN ! .Doppler width of line DOPP=WAVSY(J)*SQRT(CDOPPT*TEMP(I)/ATW(NSPSY(J))+ : CDOPPV*VT*VT) ! ,Damping parameter AV0=(ELNUM*EDAMSY(J)+RDAMSY(J)) ! .compute the line cutoff limit: sqrt(DOPCUT*FRAC*GAMMA(HWHM)) ! ! The line profile goes as gamma / delta lamda^2, so the cutoff ! in delta lambda should go as sqrt (av0). ! ! I tested various forms for the following function, and ! different values for DOPCUT, using several thousand lines ! from very weak to exceptionally strong resonance lines. ! VMAX = SQRT(FRAC*DOPCUT*AV0) ! .increment the totlat line opacity using Voigt function DO K=1,NW V0=ABS(W(K)-WAVSY(J)) IF (V0.LT.VMAX) THEN RLMU(K,I)=RLMU(K,I)+FRAC*VOIGT(AV0/DOPP,V0/DOPP)/DOPP ENDIF ENDDO ENDIF 2001 Nov 21: ATDATA Default value for gamma_nat was too large, and it was not set if the tabulated value was a dummy value between 1.E-6 and 2.38E-04 ! ! Adopt classical damping constant where none given ! According to Gray (p.207), in Delta lambda form, the ! half-half width is 0.59x1.e-4 Angstrom for all lines ! Previously we used a value four times this. ! IF(RDAM(I).LT.0.59E-04) RDAM(I) = 0.59E-04 !! 2.38E-04 2001 Nov 20: LSYNT Changed wl cutoff for metal lines from V0<50 to V0<5*AV0 ! Local constants. ! Wavelength cutoff in Doppler units for caclulating metal lines ! !! choose better algorithm for selecting full line width REAL DOPCUT PARAMETER (DOPCUT = 5.0) ! don't consider line contribution more than DOPCUT doppler widths ! from the line centre. IF (V0.LT.DOPCUT*AV0) THEN RLMU(K,I)=RLMU(K,I)+FRAC*VOIGT(AV0,V0)/DOPP ....later it may then be possible to drop the resonance line fudge... ! except for resonance line wings. ELSEIF (EXCSY(J).LT.0.01 .AND. V0.LT.5*DOPCUT) THEN RLMU(K,I)=RLMU(K,I)+FRAC*VOIGT(AV0,V0)/DOPP ENDIF 2001 Nov 15: ATSCAN Spectral resolution previously provided a minimum grid spacing of 0.01A, which is not sufficient in the UV. I have introduced a parameter "RESOLUTION" currently set to 1,000,000. The smallest wavelength interval is now given by DSMALL = WL / RESOLUTION Eventually, RESOLUTION should be put into COMMON somewhere so that it can be set by the user. 2001 Nov 20: LHI (from pld) 9. a) lhi.f modified to call double precision versions of spline and gauss-hermite wights b) nr_d_gauher.f taken from numerical recipes and put into pack/nr_d/nr_d.a - problem with convergence but probably not serious c) Test for (20,000;4.0) and vt=0, 5, 10, 20, 30 looks o.k. and is in spec_test/