SPECTRUM - bug reports and modification history 1993 - 1997. C.S.Jeffery --------------------------------------------------------------------------- Executive Summary Major changes: - U: Partition functions tidied up - HEI: improved wavelength sampling in HeI profiles (4026,4388,4471,4921) - VCS: new hydrogen data from VCS, incl. Lyman,Balmer,Paschen lines - COMMON_ALL: several arrays increased in size - S2S: a separate utility for formatting line lists from Hubeny's data - binary input data converted to ascii to enable dual O/S support - extended tabular data for diffuse HeI profiles Major bugs fixed: - ATSCAN: unit conflict resolved - OUTDIS: OPEN statement conflict resolved - LHEII: correcte partition function treatment - EWABUND: convergence failures trapped - ATDATA: changed to free format read - Kiel: argument type mismatch in SEATON Persistent problems: - CONFL: overflow in C/L ratio - LSYNT: partition functions not correct - LSYNT: continuous opacity needs improving - LSYNT: problems with HeII lines - ATSCAN: requires check for existing ATSCAN - ATCOMP: need better way of defining composition for synthesis - HEII: need ability to treat as blends with Balmer lines. - HEI: need better wavelength sampling for 4026 and 4388. ------------------------------------------------------------------------ Journal of actions 2-FEB-1993 Some perceived changes to MAIN to create SPECTRUM (a different main tool) Can these be done without disturbing MAIN ? 1. Make filename conventions standard (internally and externally) done 4-FEB-1993 and applied to all files 2. Get model filename and linelist filename internally ? 3. Disallow cycling through models ? 4. Specify options and constants more conventionally ... 2-MAR-1993 Function U (-lkiel) has been modified to use TAP_PARFN (-ltap), because it is only sensible to maintain a single version. Data for Sr, Y and Zr had been added. Function QAS (-lkiel) is redundant 4-MAR-1993 Increased paramter KI (number of atoms) to 40. Changed format statement 201 in atdata.f, 100 in outdis.f and outter.f, to allow for 10 character element names. Changed data type of ELEM from REAL*8 to CHARACTER*10. Moved ELEM from common block /IONS/ to new block /IONSC/ Shuffled order of variables in common blocks /LINE/ and /HE2/ (? forgotten its proper name) to remove misalignment warnings. Still some left ! LTE_LINES .... increased length of element identifier field from 8 to 10 characters in file atoms.d. 27-MAR-1995 ATSCAN, COMMON_ALL: Unit NSCAN added for output fropm ATSCAN (conflict with NDISC in OUTDIS) OUTDIS: OPEN statement changed to STATUS='UNKNOWN' to resolve conflict with existing logical link (SPECTRUM.TMP) defined in script Spectrum. kiel / dp: Both libraries contain u.f. dp - is p.f. routines in line kiel - link to TAP p.f. routines The first needs to be removed, then ALL tested to ensure TAP is version correct. CONFL: When calculating abundances from equivalent widtsh, the output for C/L RATIO=******* always overflows. I suspect iterative calculation of CONFL in EQUIV or similar. Seems ok in normal mode. HEI: (June 1995) The forbidden component of HeI4921 (4910) was not adequately sampled in wavelength. I increased the sampling. (Yes some stars show this component). VCS: (August 1995) Michael Lemke supplied VCS tables for MANY Lyman, Balmer and Paschen lines. I am implementing code to use these. VCS_READ, VCS2, plus a modification to HSTARK. LHI and LSYNT will also need to be modified. 3-AUG-1995 LHEII: still used constant partition functions. Changed to use variable partition functions. VCS_READ: written / tested VCS2: written / tested HSTARK: modified (eventually: VCS->BALMER, VCS2->VCS) LHI: 'lh1_new.f' written. / tested LHI used constant partition functions. Changed. LSYNT: Could leave just Balmer lines for now. LSYNT needs looking at anyway because it doesn't deal with partition functions properly and does not treat the continuous opacity properly. 7-AUG-1995 VCS_READ, VCS2, HSTARK, LHI, tested and give good results. note - Lyman-alph fails in Hdef model - Why are 'Kurucz' profiles for Halpha-delta 'better' than VCS - What do we do about natural damping and resonance broadening in hydrogen lines ? 14-AUG-1995 ATHEII added He II 4859 and 4338 to lheii_red.d 16-AUG-1995 EWABUND convergence also fails under other conditions, put a limit of 25 iterations in. 17-AUG-1995 ATDATA altered read format 101 (line data) to free format (*) in order to accommodate variable length precision in the wavelength (and other) data. Watch out for problems ... 1-DEC-1995 HEI: Changed both 4921 and 4471 wavelength grids to asymmetric (ie read in all wavelengths) so that 4910 and 4517 were both fully sampled. note - problem with extrapolation of He I profiles to high density currently switched off HeI extrapolation 12-DEC-1995 ATSCAN, LHI, HSEARCH: In order to improve spectrum synthesis I included code for dealing with VCS profiles in synthesis routines. I added some code to improve frequency sampling in sparse line regions. note - synthesis not valid for large spectral regions - can be fixed - problems including 4686 and 4861 in synthesis - problems including other HeII lines in synthesis 13-DEC-1995 SYNTH: Proposal for a new structure. subroutine SYNTH should parallel WHOLE o wavelength scale (done in ATSCAN) done o continuous opacity grid (CONT_OP) o line opacity (LIN_OP) done o radiative transfer (SYN_INTEG) o remainder as in WHOLE. ATSCAN: crashes if ATSCAN already exists - need to check ATCOMP: Proposal for new structure It should be possible to read in a file which defines the composition for a synthetic spectrum separately from having to edit the wavelength grid. ATSCAN: Proposal for modification to grid definition. Addition of wavelength points for metal lines should only proceed up to half the interval between the current line and the next nearest lines on either side. COMMON_OPAC: B star synthetic spectra were giving wrong results. Found that common_all.f had evolved beyond common_opac.f, so the qub library was entirely out of date. Easily corrected by recompiling and linking. 15-MAR-1996: S2S: I wrote a natty little utility for converting parts of Hubeny's versions of the Kurucz linelists into a format suitable for inserting into the linelist data file. It's not an ideal way to do it, but it works and allows me to do spectrum synthesis for limited sections of UV spectrum. Works very nicely, although the spectral resolution in the code is currently lower than the GHRS! KIEL opacity routines: Tried to synthesize UV spectra for Kiel model atmosphere (H-rich). I got no lines at 1335 or at 1398/1402, but did get some weak CIV lines at 1550. I got good results at 4267. Where I got no lines, TAUC had a very low value (6.3E-7 = lower limit of atmosphere) right across the line. I got good results at 1340-45 A for a QUB standard model (t25g40). I 'editted' the Kiel model to make Spectrum think it was a QUB model, and got good results at 1340-45 A (similar to above). I conclude there is a bug in the Kiel opacity routines which is setting the continuous opacity to a very high value below or around some critical frequency. By making tests with a series of dummy lines, I tracked the critical wavelength to 1420 A. The error turns out to be a call to SEATON by sb OPKCA2 where the type of constants in the argument list do not match the double precision type of the dummy arguments in SEATON. - The problem is fixed. Lyman-alpha: Still having problems with some of the new VCS profiles, including Lyman alpha. 20-MAR-1996 Lyman-alpha: LHI_ID The code for identifying lines was incorrect, and is now fixed. Lyman-alpha was still appearing much weaker than expected. LHI The code for the Boltzmann equation in LHI has the excitation potential hardwired for the Balmer series!! Changed the equation to include EH(NLOW) instead of 0.75*CHYDEV, and entered data for EH into COMMON_VCS. I have a nagging doubt about the expression for FEX, because for resonance lines the quantity exp(-CEVTOT*EXCIT(NLINE)/T) = 1, so the sum over all levels will be greater than 1... IS THIS RIGHT ? Otherwise the Ly-alpha models are now looking good. 18-APR-1997 VCS and ATHEI Three data files (balmer, prof_4471 and prof_4922) were being read as binary files. This makes portability a minor problem. The convenient solution (not the computationally most efficient) was to convert the binary profiles to ascii and to convert the code to read the ascii files. 29-APR-1997 TRANSFER, I_INTEG and F_INTEG There are problems - 1) redesign the makefile(s) so that the code builds correctly under both Solaris and Dec/Unix. 2) model : t20g40 lines : csj pld ml-vcs extrap-off 1 4101 12 5 6.262 6.262 6.260 1 4343 12 5 6.024 6.025 6.02 1 4861 12 5 5.856 5.856 5.847 1 6560 12 5 4.322 4.322 2 4471 11 5 1.641 2.738 2.738 2 4922 11 5 0.920 1.596 1.596 runs ok, although equivalent widths do not look great, they are identical on Solaris and Dec systems. My results - csj Checking equivalent widths with Philip Dufton - pld 3) model : /simon/stars/models/0099001/1725L.q lines : 1 4101 8 10 1 4343 8 10 causes run-time problems but only for Feutrier solution of transfer equation. Problem traced to alignment of variables in COMMON /source_block/ which is only referenced by subroutines TRANSFER, I_INTEG and F_INTEG. Problem was cleared. 30-APR-1997 Cleared out redundant files from all folders, especially from /data and data/utils/. The latter were put into /archive. NOTE: Spectrum only reads in a reduced set of HeII line profiles at present, including 10125, 5412, 4859, 4542, 4339, 4200 and 4686. What is it doing for blends with Hydrogen Balmer lines ? The HeII data is available in a more complete dataset. 1-MAY-1997 Numerical profiles for Balmer lines and HeI lines are wrong - see above (4343, 4722 and 4922) 2-MAY-1997 HSTARK Following a bug report from pld, corrected HSTARK to use old versions of tables for Halpha - Hepsilon, rather than Halpha - Hdelta. ATHEI Converted the disk HeI profiles to logarithmic ascii form. Saves some arithmetic, but doesn't fix the numerical problem with the equivalent widths. LHI_ID - H-gamma I changed HSTARK over so that Hbeta - Hepsilon are computed from Lemke's VCS tables - and I get exactly the same error occurring for Hgamma. Hbeta and Hdelta are OK to within errors. I had been using 4343 as the wavelength for H-gamma! When I changed it to 4340.5 (hey presto!) the answers started to come out right. On experiment I found I could put in things like 1 4097 12 5 and it would recognise it as H-delta, but give an equivalent width that was slightly incorrect. The failure comes at exactly 2A error in wavelength, but there are minor changes in equivalent width for smaller errors. I have inserted a trap in LHI_ID to detect wavelengths out by more than 1 Angstrom. A warning is given, but the program will attempt to continue without taking corrective action. This does not seem to be the cause of the problem for 4471 and 4922 5-MAY-1997 HEIDIF There is a section in HEIDIF where I switch extrapolation in electron density off in order to avoid problems (see above). The version I had running had the switching code commented out (for what reason I cannot imagine), so that I was extrapolating the HeI profiles at high electron density. When I replaced the code, and thus switched off the extrapolation I reproduced pld's results exactly. I obtained the following equivalent widths for a series of models with and without high-density extrapolation for HeI 4471. Note that WITH extrapolation, 4471 is weaker than 4026. It is important to realise that line shifts cause high-ne extrapolation to REDUCE the opacity in the line cores, which is opposite to that expected. It is not really better to avoid extrapolation; it is simply less wrong. The lines should probably be stronger than we have calculated, but we do not have the necessary theoretical profiles. I think this is IMPORTANT for B stars as well as peculiar objects like helium stars. Alain Beauchamp (1997: ApJ Supp 108, 559) has calculated profiles for 21 HeI lines and their forbidden components to high electron densities (6.E17/cm^3). The BCS tables only go to 1.E16/cm^3. However, I have not yet had success in obtaining the data from Montreal. Beachamp has left astronomy. Note that at present extrapolation is still done for 4026 and 4388 ..... Source function: Planck function. He abundance: 11.0, microturbulent velocity 5 km/s: 4471 extrap-on extrap-off 4026 4388 t20g45 2.018 3.457 2.310 1.167 t20g40 1.641 2.738 1.854 0.965 t20g35 1.267 2.078 1.388 0.751 t20g30 0.885 1.496 0.941 0.539 ...... *** 20-MAY-1997 *** I had left another bug in HEIDIF which prevented interpolation at lower densities than the maximum density tabulated. So the extrap-off version was also giving suspect results. *** EWABUND Under rare circumstances; the iteration for abundance would oscillate about a solution with a radius slightly greater than the convergence criterion. On exit from the DO loop, WHOLE would crash - for reasons not yet determined. The convergence criterion has been relaxed for the time being to give equivalent width better than 0.1% or abundance better than 1%. HEI DIFFUSE LINES - spectrum - CHECK 4026 Shamey 1969 - log ne = 14 - 16; 40,000 K 4388 Shamey 1969 - log ne = 14 - 16; 40,000 K 4471 Barnard Cooper & Smith 1974; log ne = 13 - 16; 4 temperatures OK 4921 Barnard Cooper & Smith 1975; log ne = 13 - 16; 4 temperatures - synspec - 4026 Shamey 1969 - log ne = 14 - 17.5; 4 temperatures 4388 Shamey 1969 - log ne = 14 - 17.5; 4 temperatures 4471 Barnard Cooper & Smith 1974; log ne = 13 - 16; 4 temperatures 4921 Barnard Cooper & Shamey 1969; log ne = 14 - 17.5; 4 temperatures Note: BCS69 do not print intermediate densities (15.5,16.5,... etc), but offer them for private communication. I do not have them Changes: Incorporate high density profiles for 4471 and 4921 Incorporate high densities and all temperatures for 4388 and 4026 The file names are changed as follows prof_4921.d -> he1_4921.bcs75 prof_4471.d -> he1_4471.bcs74 he1_4471.bcs74 + he1_4471.bcs69 -> prof_4471.d he1_4921.bcs75 + he1_4921.bcs69 -> prof_4921.d he1_4026.sha69 -> prof_4026.d he1_4388.sha69 -> prof_4388.d HeI line profile file formats are defined by extensions .1 and .2, thus xxx.1 - Hubeny et al. format - kappa xxx.2 - Spectrum format - log(kappa) Convert using program hei_conv.f All but the prof_files are moved to the archive folder 20-MAY-1997 COMMON_ALL, ATHEI, HEIDIF common block HEIPROF -> HEI and integrated, same procedures used for reading and interpolating in all four profiles Startlingly it works - EXCEPT 4471 and 4388 look more like 'extrapolated' results (see above). This is because I fixed a bug in the 'extrapolation' limit which was also preventing interpolation at high density. 4471 4921 4026 4388 t20g45 2.017 1.146 2.202 1.195 t20g40 1.640 0.919 1.767 0.948 t20g35 1.265 0.699 1.335 0.731 t20g30 0.883 0.499 0.915 0.520 I "think" this is now working. 21-MAY-1997 WAVEG Modified to generate better wavelength grids for 4026 and (in particular) 4388. Both have forbidden components, well-separated from the main line centre. 4383.3 shows up in computed profiles (Shamey data) for helium stars. The component at 4045.2 is not present in the Shamey tables. 3-JUNE-1997 LSYNT Modified to allow for inclusion of HeII lines. Line opacity constant changed for each line in synthesis (affected H and HeII) INTEG_FEAUT Uses local continuous opacity for spectral synthesis. It is recalculated every DELTA * lambda. The continuum flux is interpolated for every wavelength point, so that normalization is done accurately. Currently DELTA = 0.001. EQUIV New normalization ofr synthetic spectra (NMODE=10, NSOURCE=3) Other routines to change: ATSCAN - include all Balmer lines ? ATSCAN - search for HeII lines - INTEG_FEAUT gives rubbish answers if a HeII line is included in the synthesis. HESEARCH - not compatible with HEIDIF: i.e. gf-values and wavelength grids