fmq/TLUSTY
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

0

Browse Source
This commit is contained in:
fmq
2026-07-21 22:25:14 +08:00
parent 15971a5e01
commit c53fc99c1d
298 changed files with 459694 additions and 31 deletions
Download Patch File Download Diff File Expand all files Collapse all files
synspec/LINDAT.FOR
+87
View File
@@ -0,0 +1,87 @@
PARAMETER (MLIN0 =1200000,
* MGRIEM = 10,
* MNLT = 2000,
* MSPHE2 = 20,
* MLIN = 190000,
* MPRF = MLIN0)
C
PARAMETER (MLINM0 =9000000,
* MLINM =1000000,
* MMLIST = 3)
C
REAL*4 EXCL0(MLIN0),
* EXCU0(MLIN0),
* GF0(MLIN0),
* EXTIN(MLIN0),
* BNUL(MLIN0),
* GAMR0(MPRF),
* GS0(MPRF),
* GW0(MPRF),
* WGR0(4,MGRIEM),
* EXCLM(MLINM0,MMLIST),
* GFM(MLINM0,MMLIST),
* EXTINM(MLINM0,MMLIST),
* GRM(MLINM0,MMLIST),
* GSM(MLINM0,MMLIST),
* GWM(MLINM0,MMLIST),
* GVDWH2(MLINM0,MMLIST),
* GEXPH2(MLINM0,MMLIST),
* GVDWHE(MLINM0,MMLIST),
* GEXPHE(MLINM0,MMLIST)
C
COMMON/LINTOT/FREQ0(MLIN0),
* EXCL0,
* EXCU0,
* GF0,
* EXTIN,
* BNUL,
* INDAT(MLIN0),
* INDNLT(MLIN0),
* ILOWN(MLIN0),
* IUPN(MLIN0),
* IJCONT(MLIN0),
* INDLIN(MLIN),
* INDLIP(MLIN),
* NLIN0,NLIN,IRLIST,
* NNLT,NGRIEM
C
COMMON/MOLTOT/FREQM(MLINM0,MMLIST),
* EXCLM,
* GFM,
* EXTINM,
* GRM,GSM,GWM,
* GVDWH2,GEXPH2,GVDWHE,GEXPHE,
* INDATM(MLINM0,MMLIST),
* INMLIN(MLINM,MMLIST),
* INMLIP(MLINM,MMLIST),
* NLINM0(MMLIST),
* NLINML(MMLIST),
* NLINMT(MMLIST),
* IUNITM(MMLIST),
* INACTM(MMLIST),
* IVDWLI(MMLIST),
* NMLIST
CHARACTER*40 AMLIST(0:MMLIST)
COMMON/LISPAR/AMLIST,
* IBIN(0:MMLIST)
C
COMMON/LINPRF/GAMR0,
* GS0,
* GW0,
* WGR0,
* IPRF0(MPRF),
* ISPRF(MPRF),
* IGRIEM(MPRF),
* ISP0(MSPHE2),NSP
C
COMMON/LINNLT/ABCENT(MNLT,MDEPTH),
* SLIN(MNLT,MDEPTH)
C
COMMON/LINDEP/PLAN(MDEPTH),
* STIM(MDEPTH),
* EXHK(MDEPTH)
C
COMMON/LINCTR/DFRCON,IJCNTR(MLIN),IJCMTR(MLINM,MMLIST)
COMMON/MLINRE/FRLASM(MMLIST),ALASTM(MMLIST),TMLIM(MMLIST),
* NXTSEM(MMLIST),IPRSEM(MMLIST),IREADM(MMLIST)
synspec/MODELP.FOR
+65
View File
@@ -0,0 +1,65 @@
C
C Basic parameters of the model atmosphere
C
COMMON/MODELP/DM(MDEPTH),
* TEMP(MDEPTH),
* ELEC(MDEPTH),
* DENS(MDEPTH),
* ZD(MDEPTH),
* VTURB(MDEPTH),VTB,
* ABSTD(MDEPTH),
* ABSTDW(MFREQC,MDEPTH),
* POPUL(MLEVEL,MDEPTH),
* POPREL(MLEVEL,MDEPTH),
* DMR0(MDEPTH),
* DMRP(MDEPTH),
* SBF(MLEVEL),
* USUM(MIOEX),
* WOP(MLEVEL,MDEPTH),
* WNHINT(NLMX,MDEPTH),
* WNHE2(NLMX,MDEPTH),
* RRR(MDEPTH,MION,MATOM),
* JT(MDEPTH),
* TI0(MDEPTH),
* TI1(MDEPTH),
* TI2(MDEPTH)
character*8 cmol(mmolec)
COMMON/MOLPAR/RRMOL(MMOLEC,MDEPTH),
* DOPMOL(MMOLEC,MDEPTH),
* AMMOL(MMOLEC),
* CMOL,
* anh2(mdepth),anch(mdepth),anoh(mdepth),
* anhm(mdepth)
C
COMMON/OPACAT/OPATM(MATOM,MFREQ,MDEPTH),
* EMATM(MATOM,MFREQ,MDEPTH),
* OPATML(MATOM,MFREQ),
* GRADAT(MATOM,MDEPTH),
* GRADFA(MATOM,MDEPTH),
* POPAT(MATOM,MDEPTH),
* DGRAD0(MATOM,MATOM,MDEPTH),
* DGRADP(MATOM,MATOM,MDEPTH)
C
COMMON/RADFLD/RAD(MFREQ,MDEPTH),
c * FAK(MFREQ,MDEPTH),
c * ALI(MFREQ,MDEPTH),
c * FLXH(MFREQ,MDEPTH),
* RAD0(MFREQ,MDEPTH),
* FLX0(MFREQ,MDEPTH),
* flxt(mdepth),
* flxi(mdepth)
C
COMMON/XENPRF/PRFXB(MLINH,MHWL,MHT,MHE),
* PRFXR(MLINH,MHWL,MHT,MHE),
* PRFB(MLINH,MDEPTH,MHWL),
* PRFR(MLINH,MDEPTH,MHWL),
* ALXEN(MLINH,MHWL),
* XTXEN(MHT,MLINH),
* XNEXEN(MHE,MLINH),XNEMIN,
* NWLXEN(MLINH),
* NTHXEN(MLINH),
* NEHXEN(MLINH),
* ILXEN(4,22),
* IHXENB
C
synspec/OPTPAR.FOR
+4
View File
@@ -0,0 +1,4 @@
PARAMETER (MFRTAB = 100000,
* MTTAB = 20,
* MRTAB = 20,
* MSFTAB = 2000000.
synspec/PARAMS.FOR
+223
View File
@@ -0,0 +1,223 @@
C
C Parameters that specify dimensions of arrays
C
IMPLICIT REAL*8 (A-H, O-Z),LOGICAL*1 (L)
character*4 typat
PARAMETER (MATEX = 30,
* MIOEX = 90,
* MLEVEL= 1650,
* MDEPTH= 100,
* MDEPF = 500,
* MFREQ = 2000,
c * MFREQ = 120,
* MFREQC= 2000,
* MFRQ = 2000,
* MOPAC = MFRQ,
* MMU = 20,
* MCROSS= MLEVEL,
* MFIT = 1650,
* MFCRA = 1200,
* MTRAD = 3,
* MATOM = 99,
* MATOMBIG = 99,
* MION = 90,
* MION0 = 9,
* MMOLEC=500,
* MPHOT = 10,
* MZZ = 2,
* MMER = 2,
* NLMX = 80,
* MI1 = MION0-1,
* MLINH = 78,
* MHT = 7,
* MHE = 20,
* MHWL = 55)
PARAMETER (MFGRID = 100000,
* MTTAB = 21,
* MRTAB = 20,
* MSFTAB = 6000000)
parameter (mfhtab=1000,
* mtabth=10,
* mtabeh=10)
c
C Basic physical constants
C
PARAMETER (H = 6.6256D-27,
* CL = 2.997925D10,
* BOLK = 1.38054D-16,
* HK = 4.79928144D-11,
* EH = 2.17853041D-11,
* BN = 1.4743D-2,
* SIGE = 6.6516D-25,
* PI4H = 1.8966D27,
* HMASS = 1.67333D-24)
C
C Unit number
C
PARAMETER (IBUFF=95)
C
C Variables to hold quantum numbers limits
C (see LEVLIMITS below)
C
INTEGER*4 SQUANT1(MLEVEL),SQUANT2(MLEVEL),
* LQUANT1(MLEVEL),LQUANT2(MLEVEL),
* PQUANT1(MLEVEL),PQUANT2(MLEVEL)
C
C Basic parameters
C
COMMON/BASNUM/NATOM,
* NION,
* NLEVEL,
* ND,NDSTEP,
* NFREQ,NFROBS,NFREQC,NFREQS,
* NMU
COMMON/LTESET/LTE,LTEGR
COMMON/INPPAR/TEFF,
* GRAV,
* YTOT(MDEPTH),
* WMM(MDEPTH),
* WMY(MDEPTH),
* vaclim,
* ATTOT(MATOM,MDEPTH)
COMMON/BASICM/IMODE,
* IMODE0,
* IFREQ,
* INLTE,
* IDSTD,
* IFWIN,
* IFEOS,
* IBFAC
COMMON/INTKEY/INMOD,INTRPL,ICHANG,ICHEMC,IATREF,ICONTL
COMMON/LBLANK/IBLANK,NBLANK
COMMON/NXTINI/ALM00,ALST00,NXTSET,INLIST,ALAMBE,DLAMLO
COMMON/IPRNTR/IPRIN
C
C Parameters for explicit atoms
C
COMMON/ATOPAR/AMASS(MATEX),
* ABUND(MATEX,MDEPTH),
* RELAB(MATEX,MDEPTH),
* NUMAT(MATEX),
* N0A(MATEX),
* NKA(MATEX),
* SABND(MATEX)
C
C Parameters for explicit ions
C
COMMON/IONPAR/FF(MIOEX),
* NFIRST(MIOEX),
* NLAST(MIOEX),
* NNEXT(MIOEX),
* IUPSUM(MIOEX),
* IZ(MIOEX),
* IFREE(MIOEX),
* INBFCS(MIOEX),
* ILIMITS(MIOEX)
C
C Parameters for explicit levels
C
COMMON/LEVPAR/ENION(MLEVEL),
* G(MLEVEL),
* NQUANT(MLEVEL),
* IATM(MLEVEL),
* IEL(MLEVEL),
* ILK(MLEVEL),
* ifwop(mlevel),
* isemex(matom)
C
C Limits for explicit levels
C
COMMON/LEVLIMITS/ENION1(MLEVEL),
* ENION2(MLEVEL),
* SQUANT1,
* SQUANT2,
* LQUANT1,
* LQUANT2,
* PQUANT1,
* PQUANT2
C
C Parameters for all considered transitions
C
COMMON/TRAPAR/IBF(MLEVEL),
* S0BF(MLEVEL),
* ALFBF(MLEVEL),
* BETBF(MLEVEL),
* GAMBF(MLEVEL)
C
COMMON/MRGPAR/SGM0(MMER),
* FRCH(MMER),
* SGEXT1(MMER,MDEPTH),
* GMER(MMER,MDEPTH),
* SGMSUM(NLMX,MMER,MDEPTH),
* SGMG(MMER,MDEPTH),
* IMRG(MLEVEL),
* IIMER(MMER)
C
COMMON/DWNPAR/ELEC23(MDEPTH),
* Z3(MZZ),
* DWC1(MZZ,MDEPTH),
* DWC2(MDEPTH)
C
C additional opacities
c
COMMON/OPCPAR/IOPADD,
* IOPHMI,
* IOPH2P,
* IOPHEM,
* IOPCH,
* IOPOH,
* IOPH2M,
* IOH2H2,IOH2HE,IOH2H1,IOHHE,
* IOPHLI,
* IRSCT,
* IRSCHE,
* IRSCH2
C
C Auxiliary parameters
C
COMMON/AUXIND/IATH,IELH,IELHM,N0H,N1H,NKH,N0HN,N0M,
* IATHE,IELHE1,IELHE2
COMMON/MOLFLG/TMOLIM,MOLIND(11000),NMOLEC,IFMOL,
* MOLTAB,IRWTAB,IIRWIN,IPFEXO
COMMON/QFLAGS/ERANGE,ISPICK,ILPICK,IPPICK
C
C Parameters for atoms considered in line blanketing opacity
C
LOGICAL LGR(MATOM),LRM(MATOM)
COMMON/PFSTDS/PFSTD(MION,MATOM),MODPF(MATOM)
COMMON/ADDPOP/RR(MATOM,MION)
COMMON/ATOBLN/ENEV(MATOM,MI1),AMAS(MATOM),ABND(MATOM),
* ABNDD(MATOM,MDEPTH),ABNREF(MDEPTH),TYPAT(MATOM),
* IATEX(MATOM),INPOT(MATOM,MION0)
COMMON/ATOINI/NATOMS,IONIZ(MATOM),LGR,LRM
c
c parameters for hydrogen Stark broadening tables
c
COMMON/HYDPRF/PRFHYD(MLINH,MDEPTH,MHWL),
* WLHYD(MLINH,MHWL),
* NWLHYD(MLINH),
* WL(MHWL,MLINH),
* XT(MHT,MLINH),
* XNE(MHE,MLINH),
* PRF(MHWL,MHT,MHE,MLINH),
* WLINE(4,22),
* NWLH(MLINH),
* NTH(MLINH),
* NEH(MLINH),
* ILIN0(4,22),
* ILEMKE,
* NLIHYD
COMMON/AUXHYD/XK,FXK,BETAD,DBETA,BERGFC,CUTLYM,CUTBAL
COMMON/HHEPRF/IHYDPR,IHE1PR,IHE2PR
COMMON/HYLPAR/IHYL,ILOWH,M10,M20
COMMON/HYLPAW/IHYLW(MFREQ),ILOWHW(MFREQ),
* M10W(MFREQ),M20W(MFREQ)
COMMON/HE2PAR/IFHE2,IHE2L,ILWHE2,MHE10,MHE20
COMMON/HE2PAW/IHE2LW(MFREQ),ILWHEW(MFREQ),
* MHE10W(MFREQ),MHE20W(MFREQ)
C
C parameters for the macroscopic velocity field and angles
C
COMMON/VELPAR/ANGL(MMU),WANGL(MMU),VELC(MDEPTH),NMU0,IFLUX
synspec/SYNTHP.FOR
+10
View File
@@ -0,0 +1,10 @@
COMMON/FREQSY/FREQ(MFREQ),W(MFREQ),WLAM(MFREQ),
* FRX1(MFREQ),FRX2(MFREQ),BNUE(MFREQ),
* FRQOBS(MFREQ),WLOBS(MFREQ),
* FREQC(MFREQC),WLAMC(MFREQC),
* IJCINT(MFREQ)
COMMON/CRSAVG/FRECR(MCROSS,MFCRA),CROSR(MCROSS,MFCRA),
* CRMX(MCROSS),NFCR(MCROSS),IASV
COMMON/CRSAVQ/FRECQ(MPHOT,MFCRA),QHOT(MPHOT,MFCRA),
* AQHT(MPHOT),EQHT(MPHOT),GQHT(MPHOT),
* CRMY(MPHOT),NFQHT(MPHOT),NQHT
synspec/WINCOM.FOR
+15
View File
@@ -0,0 +1,15 @@
PARAMETER (MRCORE=20,
* MKU=MDEPTH+MRCORE,
* MEXT=MKU)
COMMON/COMANG/BMU(MKU,MDEPTH),WMUJ(MKU,MDEPTH),WMUH(MKU)
COMMON/CORADI/RD(MDEPTH),RCORE,RFNORM,PIM(MKU),RAD1(MDEPTH),
* DELZ(MKU,MDEPTH),NUD(MKU),NUDF(MKU),KMU,NREXT,
* NRCORE,NFIRY,NDF
COMMON/CORAF/DELZF(MEXT,MDEPF ),DFRQF(MEXT,2*MDEPF )
COMMON/COVEL/VEL(MDEPTH),DFRQ(MKU,2*MDEPTH),DVD(MDEPTH),
* XMDOT,XMD4,BETAV,VINF
COMMON/EXTMOD/FFQ(MOPAC),FFQV(MOPAC),RDF(MDEPF ),DENSF(MDEPF ),
* VELF(MEXT,MDEPF ),DRAY(MEXT,2*MDEPF ),
* KRAY(MEXT,2*MDEPF ),NOPAC
COMMON/OPAVEL/WDIL(MDEPTH),PLANW(MDEPTH),TRAD(MTRAD,MDEPTH),
* DENSCON(MDEPTH)
synspec/ngrid.f
+288
View File
@@ -0,0 +1,288 @@
program ngrid
C =============
c
IMPLICIT REAL*8 (A-H, O-Z)
PARAMETER (MFGRID = 300000,
* MTTAB = 25,
* MRTAB = 15,
* MSFTAB = 6000000)
c
real*4 absgrd(mttab,mrtab,mfgrid)
common/gridp0/tempg(mttab),densg(mrtab),elecgr(mttab,mrtab),
* densg0(mttab),temp1,ntemp,ndens
common/gridf0/wlgrid(mfgrid),nfgrid
common/fintab/absgrd
common/initab/absop(msftab),wltab(msftab)
common/elecm0/elecm(mttab)
dimension tempvec(mttab),rhovec(mrtab)
dimension abgrd(mfgrid),abunt(30),abuno(30)
dimension typa(30)
character*(80) tabname,optable
character*4 typa
c
istop=0
read(5,*) nlamb,inttab,wlam1,wlam2
read(5,*) tabname,ibingr
if(nlamb.gt.mfgrid) then
write(*,*) 'nlambda.gt.mfgrid - recompile with larger mfgrid'
istop=1
end if
ibinop=0
read(5,*,err=5,end=5) optable,ibinop
5 continue
write(*,*) 'original optable ',optable
c
nfgrid=nlamb
wl1=log(wlam1)
wl2=log(wlam2)
dwl=(wl2-wl1)/(nfgrid-1)
do i=1,nfgrid
wlgrid(i)=exp(wl1+(i-1)*dwl)
end do
c
c read the header of the old opacity table
c
if(ibinop.eq.0) then
open(52,file=optable,status='old')
else
open(52,file=optable,form='unformatted',status='old')
end if
c
if(ibinop.eq.0) then
read(52,*)
read(52,*)
do iat=1,30
read(52,*) typa(iat),abunt(iat),abuno(iat)
end do
read(52,*)
read(52,*)
read(52,*) ifmolt,tmolit
read(52,*)
read(52,*)
read(52,*) iophmt,ioph2t,iophet,iopcht,iopoht
read(52,*)
read(52,*)
read(52,*) numfre0,numtem0,numrh0
read(52,*)
read(52,*) (tempvec(i),i=1,numtem0)
read(52,*)
read(52,*) (rhovec(j),j=1,numrh0)
read(52,*)
read(52,*) ((elecgr(i,j),j=1,numrh0),i=1,numtem0)
else
do iat=1,30
read(52) typa(iat),abunt(iat),abuno(iat)
end do
read(52) ifmolt,tmolit
read(52) iophmt,ioph2t,iophet,iopcht,iopoht
read(52) numfre0,numtem0,numrh0
read(52) (tempvec(i),i=1,numtem0)
read(52) (rhovec(j),j=1,numrh0)
read(52) ((elecgr(i,j),j=1,numrh0),i=1,numtem0)
end if
ntemp=numtem0
ndens=numrh0
write(*,*) 'ntemp,ndens',ntemp,ndens
if(ntemp.gt.mttab) then
write(*,*) 'ntemp.gt.mttab - recompile with larger mttab'
istop=1
end if
if(ndens.gt.mrtab) then
write(*,*) 'ndens.gt.mrtab - recompile with larger mrtab'
istop=1
end if
c
if(istop.eq.1) stop
c
c read file fort.27 (internal opacitis) and set the
c arrays for the new opacity tables
c
if(inttab.ne.1) then
c
c 1st possibility - approximate preserving an integral of opacity
c
indext=0
indexn=0
nfr=0
10 continue
indext=indext+1
20 continue
indexn=indexn+1
nfr=0
ijgrd=0
30 continue
ijgrd=ijgrd+1
wlgr=0.5*(wlgrid(ijgrd)+wlgrid(ijgrd+1))
isum=0
sum=0.
40 continue
read(27,*,err=50,end=50) ip,wl,abl
if(wl.le.wlgr) then
sum=sum+exp(abl)
isum=isum+1
if(ip.le.nfr) go to 50
nfr=ip
go to 40
end if
abgrd(ijgrd)=log(sum/float(isum))
if(ijgrd.lt.nfgrid) go to 30
50 continue
write(*,*) 'it,ir,nf',indext,indexn,nfr,nfgrid
do ij=1,nfgrid
absgrd(indext,indexn,ij)=real(abgrd(ij))
end do
absgrd(indext,indexn,nfgrid)=absgrd(indext,indexn,nfgrid-1)
c
if(indexn.lt.ndens) then
go to 20
else
indexn=0
end if
if(indext.lt.ntemp) then
go to 10
else
indext=0
end if
c
else
c
c 2nd possibility - an interpolation in wavelengths
c
indext=0
indexn=0
nfr=0
110 continue
indext=indext+1
120 continue
indexn=indexn+1
nfr=0
130 continue
read(27,*,err=150,end=150) ip,wl,abl
if(ip.le.nfr) go to 150
wltab(ip)=wl
absop(ip)=abl
nfr=ip
go to 130
150 continue
write(*,*) 'it,ir,nf',indext,indexn,nfr,nfgrid
call intrp(wltab,absop,wlgrid,abgrd,nfr,nfgrid)
do ij=1,nfgrid
absgrd(indext,indexn,ij)=real(abgrd(ij))
end do
c
if(indexn.lt.ndens) then
go to 120
else
indexn=0
end if
if(indext.lt.ntemp) then
go to 110
else
indext=0
end if
end if
c
c store the opacities in the new table
c
if(ibingr.eq.0) then
open(53,file=tabname,status='unknown')
write(53,600)
do iat=1,30
write(53,601) typa(iat),abunt(iat),abuno(iat)
end do
write(53,602) ifmolt,tmolit
write(53,603) iophmt,ioph2t,iophet,iopcht,iopoht
write(53,611) nfgrid,ntemp,ndens
write(53,612) (tempvec(i),i=1,ntemp)
write(53,613) (rhovec(j),j=1,ndens)
write(53,614) ((elecgr(i,j),j=1,ndens),i=1,ntemp)
do k = 1, nfgrid
write(53,615) k,wlgrid(k),2.997925e18/wlgrid(k)
do j = 1,ndens
write(53,616) (absgrd(i,j,k),i=1,ntemp)
end do
end do
600 format('opacity table with element abundances:'/
* 'element for EOS for opacities')
601 format(' ',a4,1p2e12.3)
602 format(/'molecules - ifmol,tmolim:'/,i4,f10.1)
603 format('additional opacities'/' H- H2+ He- CH OH'/
* 5i4)
611 format(/'number of frequencies, temperatures, densities:'
* /10x,3i10)
612 format('log temperatures'/(6F11.6))
613 format('log densities'/(6F11.6))
614 format('log electron densities from EOS'/(6f11.6))
615 format(/' *** frequency # : ',i8,f15.5/1pe20.8)
616 format((1p6e14.6))
else
do iat=1,30
write(63) typa(iat),abunt(iat),abuno(iat)
end do
write(63) ifmolt,tmolit
write(63) iophmt,ioph2t,iophet,iopcht,iopoht
write(63) nfgrid,ntemp,ndens
write(63) (tempvec(i),i=1,ntemp)
write(63) (rhovec(j),j=1,ndens)
write(63) ((elecgr(i,j),j=1,ndens),i=1,ntemp)
do k = 1, nfgrid
write(63) wlgrid(k)
do j = 1, ndens
write(63) (absgrd(i,j,k),i=1,ntemp)
end do
end do
end if
c
end
C
C
C ****************************************************************
C
C
subroutine intrp(wltab,absop,wlgrid,abgrd,nfr,nfgrid)
c =====================================================
c
IMPLICIT REAL*8 (A-H, O-Z)
PARAMETER (MFGRID = 400000)
c
dimension wltab(1),absop(1),wlgrid(1),abgrd(1)
dimension yint(mfgrid),jint(mfgrid)
c
c set up interpolation coefficients for frequency interpolation
c by bisection
c
fr1=wltab(1)
fr2=wltab(nfr)
do ij=1,nfgrid
xint=wlgrid(ij)
jl=0
ju=nfr+1
10 continue
if(ju-jl.gt.1) then
jm=(ju+jl)/2
if((fr2.gt.fr1).eqv.(xint.gt.wltab(jm))) then
jl=jm
else
ju=jm
end if
go to 10
end if
j=jl
if(j.eq.nfr) j=j-1
if(j.eq.0) j=j+1
jint(ij)=j
c yint(ij)=un/log10(wltab(j+1)/wltab(j))
yint(ij)=1./(wltab(j+1)-wltab(j))
end do
c
do ij=1,nfgrid
j=jint(ij)
rc=(absop(j+1)-absop(j))*yint(ij)
c abgrd(ij)=rc*log10(wlgrid(ij)/wltab(j))+absop(j)
abgrd(ij)=rc*(wlgrid(ij)-wltab(j))+absop(j)
end do
c
return
end
synspec/rotin.f
+730
View File
@@ -0,0 +1,730 @@
PROGRAM ROTIN3
c
c -------------------------------------------------------------------
c Program for performing rotational and instrumental convolution
c for a calculated stellar spectrum obtained by program SYNSPEC
c -------------------------------------------------------------------
c
PARAMETER (MLAM=1000000,
* MCON=20000)
PARAMETER (MCONV=50000)
implicit real*8 (a-h,o-z)
DIMENSION FLAM0(MLAM),FLAM1(MLAM),FLAM2(MLAM),FCN(MLAM),
* WLAM0(MLAM),WLAM1(MLAM),WLAM2(MLAM),
* FCON0(MCON),WCON0(MCON),
* wlam3(mlam),flam3(mlam)
character*20 fname7,fname17,fnout
c EQUIVALENCE (FLAM2(1),FLAM0(1)),
c * (WLAM2(1),WLAM0(1)),
c * (FCN(1),FLAM1(1))
T0=0.
c
c ---------------------------------------------------------------
C INPUT - from unit 5 - four lines of input
c ---------------------------------------------------------------
c
c 1. filenames:
c
c fname7 - name of the file containing the detailed synthetic spectrum
c (i.e. fort.7 produced by Synspec )
c
c fname17 - name of the file containing the continuum flux
c (i.e. fort.17 produced by Synspec )
c
c fnout - name of the output file - convolved spectrum
c
c
c 2. parameters for rotational convolution
c
c VROT - v sin i (in km/s)
c if VROT=0 - rotational convolution is
c a) either not calculated,
c b) or, if simultaneously FWHM is rather large
c (vrot/c*lambda < FWHM/20.),
c vrot is set to FWHM/20*c/lambda;
c if VROT >0 but the previous condition b) applies, the
c value of VROT is changed as in the previous case
c if VROT<0 - the value of abs(VROT) is used regardless of
c how small compared to FWHM it is
c CHARD - characteristic scale of the variations of unconvolved
c stellar spectrum (basically, characteristic distance
c between two neighbouring wavelength points) - in A
c - if =0 - program sets up default (0.01 A)
c STEPR - wavelength step for evaluation rotational convolution;
c - if =0, the program sets up default (the wavelength
c interval corresponding to the rotational velocity
c devided by 3.)
c if <0, convolved spectrum calculated on the original
c (detailed) SYNSPEC wavelength mesh
c
c
c 3. parameters for instrumental convolution
c
c FWHM - full width at half maximum for Gaussian instrumental
c profile
c STEPI - wavelength step for evaluating instrumental convolution
c - if =0, the program sets up default (FWHM/10.)
c - if <0, convolved spectrum calculated with the previous
c wavelength mesh:
c either the original (SYNSPEC) one if vrot=0,
c or the one used in rotational convolution (vrot > 0)
c
c
c 4. wavelength interval and normalization of spectra
c
c ALAM0 - initial wavelength
c ALAM1 - final wavelength
c IREL - for =1 relative spectrum
c =0 absolute spectrum
c
read(5,*) fname7,fname17,fnout
read(5,*) vrot,chard,stepr
read(5,*) fwhm,stepi,vmac
read(5,*) alam0,alam1,irel
c
open(unit=37,file=fname7,status='old')
open(unit=47,file=fname17,status='old')
open(unit=11,file=fnout)
c
c if some of the above convolution are zero, set up defaults
c
SLAM=0.5*(ALAM0+ALAM1)
DROT=SLAM*ABS(VROT)/3.E5
dins=fwhm/20.
if(chard.le.0.) chard=0.01
if(drot.lt.dins.and.vrot.ge.0.) then
drot=dins
vrot=drot*3.e5/slam
end if
if(vrot.lt.0.) vrot=abs(vrot)
if(stepr.eq.0.) stepr=max(drot/5.,chard)
if(stepi.eq.0.) stepi=max(fwhm/10.,chard)
if(stepi.lt.stepr) stepi=stepr
write(6,601) vrot,chard,stepr,fwhm,stepi,alam0,alam1,irel
601 format(' vrot =',f10.3/
* ' chard =',f10.3/
* ' step(rot.conv) =',f10.3/
* ' fwhm =',f10.3/
* ' step(inst.conv) =',f10.3/
* ' lambda initial =',f10.3/
* ' lambda final =',f10.3/
* ' irel =',i6/)
c
c -------------------------------------------------------------------
c read in and transform the output from SYNSPEC
c -------------------------------------------------------------------
c
c WLAM0 - array of wavelengths
c FLAM0 - array of fluxes
c WCON0 - arrray of wavelengths for continuum
c FCON0 - array of continuum fluxes
c
I=0
ILAM=0
ICON=0
ALM=0.
10 READ(37,*,END=20) AL,FL
IF(AL.LT.ALAM0) GO TO 10
IF(AL.GT.ALAM1) GO TO 20
ILAM=ILAM+1
IF(ILAM.GT.MLAM) GO TO 20
WLAM0(ILAM)=AL
FLAM0(ILAM)=FL
GO TO 10
20 CONTINUE
30 READ(47,*,END=40) AL,FL
IF(AL.LT.ALAM0-10.) GO TO 30
IF(AL.GT.ALAM1+10.) GO TO 40
IF(AL.LE.ALM) GO TO 30
ICON=ICON+1
IF(ICON.GT.MCON) GO TO 40
WCON0(ICON)=AL
FCON0(ICON)=FL
ALM=AL
GO TO 30
40 CONTINUE
NLAM0=ILAM
NCON0=ICON
IF(NLAM0.GT.MLAM) NLAM0=MLAM
IF(NCON0.GT.MCON) NCON0=MCON
IF(NLAM0.LE.0) GO TO 1000
IF(NCON0.LE.0.AND.IREL.GT.0) GO TO 1000
c
c if reguired, normalization to get a relative spectrum
c
if(irel.eq.1) then
call interp(wcon0,fcon0,wlam0,fcn,ncon0,nlam0,2,0,0)
do 50 i=1,nlam0
50 flam0(i)=flam0(i)/fcn(i)
end if
c
c -----------------------------------------------------------------
c Rotational convolution
c (no rotational convolution if vrot=0)
c -----------------------------------------------------------------
c
if(vrot.gt.0.) then
SLAM=0.5*(WLAM0(1)+WLAM0(NLAM0))
XN=SLAM*VROT/3.E5/CHARD
NROT=INT(XN)
IF(NROT.LT.10) NROT=10
IF(NROT.GT.MCONV) then
write(6,*) 'nrot too large'
stop
end if
if(stepr.le.0.) then
nlam1=nlam0
do 80 i=1,nlam1
80 wlam1(i)=wlam0(i)
else
XN=(ALAM1-ALAM0)/STEPR
NLAM1=INT(XN)+1
DO 90 I=1,NLAM1
90 WLAM1(I)=ALAM0+STEPR*(I-1)
end if
CALL ROTINS(1,FLAM0,FLAM1,WLAM0,WLAM1,NLAM0,NLAM1,
* NROT,VROT,0.d0)
c
else
c
c no rotational convolution if VROT=0
c
nlam1=nlam0
do 100 i=1,nlam1
wlam1(i)=wlam0(i)
flam1(i)=flam0(i)
100 continue
end if
c
c ------------------------------------------------------------------
c anisotropic macroturbulent velocity (if vmac gt 0)
c after Grey - radial-tangential macroturbulence
c -----------------------------------------------------------------
c
if(vmac.gt.0.) then
nlam3=nlam1
do i=1,nlam3
wlam3(i)=wlam1(i)
end do
nins=20
call rotins(3,flam1,flam3,wlam1,wlam3,nlam1,nlam3,
* nins,vmac,0.d0)
do i=1,nlam3
flam1(i)=flam3(i)
end do
end if
c ------------------------------------------------------------------
c Instrumental convolution
c no instrumental convolution for FWHM=0
c ------------------------------------------------------------------
c
if(fwhm.gt.0.) then
IF(VROT.LE.0.) THEN
XN=6.*FWHM/CHARD
ELSE
XN=MIN(1.8E6*FWHM/VROT/SLAM,6.*FWHM/CHARD)
END IF
NINS=INT(XN)
IF(NINS.LT.10) NINS=10
c write(6,*) 'number of integration points - instrum.',nins
IF(NINS.GT.MCONV) then
write(6,*) 'nins too large'
stop
end if
if(stepi.le.0.) then
nlam2=nlam1
do 110 i=1,nlam2
110 wlam2(i)=wlam1(i)
else
XNI=(ALAM1-ALAM0)/STEPI
NLAM2=INT(XNI)+1
DO 120 I=1,NLAM2
120 WLAM2(I)=ALAM0+STEPI*(I-1)
end if
CALL ROTINS(2,FLAM1,FLAM2,WLAM1,WLAM2,NLAM1,NLAM2,
* NINS,0.d0,FWHM)
else
nlam2=nlam1
do i=1,nlam2
wlam2(i)=wlam1(i)
flam2(i)=flam1(i)
end do
end if
c
c -----------------------------------------------------------------
c
c output of final, convolved, synthetic spectrum on standard output
c in a simple format:
c
c values of wavelengths (in A) - WLAM2 versus
c corresponding values of fluxes - FLAM2
c
do 200 i=1,nlam2
write(11,600) wlam2(i),flam2(i)
200 continue
600 format(f10.3,1pe10.3)
c
c rewind 7
c rewind 17
c
1000 continue
c
stop
end
c
c ********************************************************************
c
SUBROUTINE ROTINS(MODE,HINP,HOUT,XLAM,YLAM,NLAMX,NLAMY,
* NR,VROT,FWHM)
C
C ---------------------------------------------------------------
C
C Rotational and/or instrumental convolution
C
C MODE - for MODE = 1 - rotational convolution;
C for MODE = 2 - instrumental convolution,
C with a Gaussian instrumental profile
C HINP - array input flux (arbitrary units)
C XLAM - array of input wavelengths (in angstroms)
C HOUT - array of output (convolved) flux
C YLAM - array of wavelengths in which the output flux is calculated
C NLAMX - number of input wavelengths
C NLAMY - number of output wavelengths
C NR - number of integration points for evaluating the
C convolution integrals
C VROT - v sin i [in km/s]
C FWHM - full width at half maximum of the instrum. profile function
C ---------------------------------------------------------------------
C
implicit real*8 (a-h,o-z)
PARAMETER (MCONV=5000,
* MCONV1=MCONV+1)
PARAMETER (MLAM=100000)
PARAMETER (ONE=1.,
* TWO=2.,
* HALF=0.5)
PARAMETER (DLROT=10.)
DIMENSION HINP(1),HOUT(1),XLAM(1),YLAM(1),G(MCONV1),
* IGCALC(MLAM)
C
NR1=NR+1
DO 10 I=1,NLAMY
IGCALC(I)=0
10 CONTINUE
C
IF(MODE.EQ.1.or.mode.eq.3) THEN
DLAM=YLAM(NLAMY)-YLAM(1)
IF(DLAM.LE.DLROT) THEN
SLAM=HALF*(YLAM(1)+YLAM(NLAMY))
ELSE
NCALG=INT(DLAM/DLROT)+1
NSTEP=NLAMY/NCALG
DO 20 I=1,NCALG
IGCALC(I*NSTEP+1)=1
20 CONTINUE
SLAM=YLAM(1)
END IF
ELSE IF(MODE.EQ.2) THEN
END IF
c
c initial kernel function (rotation);
c or the general kernel function (for instrumental)
c
CALL KERNEL(MODE,NR,VROT,SLAM,FWHM,XLMAX,G)
DLAM=XLMAX/NR
c
c determine the exterior intervals
c a) beginning of the interval
c
INTR0=0
IEND0=0
X0=XLAM(1)
X1=X0+XLMAX
HM0=HINP(1)
DO 50 I=1,NLAMY
IF(YLAM(I).LE.X0) THEN
IEND0=I
HOUT(I)=HM0
ELSE IF(YLAM(I).LE.X1) THEN
INTR0=I
ELSE
GO TO 60
END IF
50 CONTINUE
60 CONTINUE
c
c b) end of the interval
c
INTR1=NLAMY
IEND1=NLAMY
X0=XLAM(NLAMX)
IF(MODE.EQ.1.AND.YLAM(NLAMY)-YLAM(1).GT.DLROT) THEN
XLMAX=YLAM(NLAMY)*VROT/2.997925E5
END IF
X1=X0-XLMAX
HP0=HINP(NLAMX)
DO 70 I=NLAMY,1,-1
IF(YLAM(I).GE.X0) THEN
IEND1=I
HOUT(I)=HP0
ELSE IF(YLAM(I).GE.X1) THEN
INTR1=I
ELSE
GO TO 80
END IF
70 CONTINUE
80 CONTINUE
C
C ------------------------------------------------------------
C wavelength by wavelength convolution; integral calculated by
C the trapezoidal rule
C ------------------------------------------------------------
C
C 1. points near the beginning of the interval
C
IF(INTR0.GE.1) THEN
K0=1
DO 390 I=IEND0+1,INTR0
HOUT(I)=0.
DO 310 K=K0,NLAMX
K0=K
IF(XLAM(K).GT.YLAM(I)) GO TO 320
310 CONTINUE
320 K0=K0-1
DO 350 J=1,NR1
A2=(J-1)*DLAM
ALAM=YLAM(I)+A2
K=K0+1
330 CONTINUE
IF(ALAM.LT.XLAM(K)) THEN
HPLUS=HINP(K-1)+(HINP(K)-HINP(K-1))/(XLAM(K)-XLAM(K-1))*
* (ALAM-XLAM(K-1))
HOUT(I)=HOUT(I)+HPLUS*G(J)
ELSE IF(ALAM.EQ.XLAM(K)) THEN
HOUT(I)=HOUT(I)+HINP(K)*G(J)
ELSE
K=K+1
GO TO 330
END IF
350 CONTINUE
C
DO 380 J=1,NR1
A2=(J-1)*DLAM
ALAM=YLAM(I)-A2
K=K0
360 CONTINUE
IF(ALAM.GT.XLAM(K)) THEN
HMINUS=HINP(K)+(HINP(K+1)-HINP(K))/(XLAM(K+1)-XLAM(K))*
* (ALAM-XLAM(K))
HOUT(I)=HOUT(I)+HMINUS*G(J)
ELSE IF(ALAM.EQ.XLAM(K)) THEN
HOUT(I)=HOUT(I)+HINP(K)*G(J)
ELSE
K=K-1
IF(K.LT.1) THEN
HOUT(I)=HOUT(I)+HM0*G(J)
GO TO 380
END IF
GO TO 360
END IF
380 CONTINUE
IF(K0.LE.0) K0=1
390 CONTINUE
END IF
C
C 2. inner points
C
c 601 format(5i5)
if(intr0.le.0) intr0=1
K0=1
DO 300 I=INTR0+1,INTR1-1
C
C re-evaluate the kernel function if necessary
c
IF(IGCALC(I).EQ.1) THEN
SLAM=YLAM(I)
CALL KERNEL(MODE,NR,VROT,SLAM,FWHM,XLMAX,G)
DLAM=XLMAX/NR
END IF
c
c perform the convolution integral
c
HOUT(I)=0.
DO 110 K=K0,NLAMX
K0=K
IF(XLAM(K).GT.YLAM(I)) GO TO 120
110 CONTINUE
120 K0=K0-1
DO 200 J=1,NR1
A2=(J-1)*DLAM
ALAM=YLAM(I)+A2
K=K0+1
130 CONTINUE
IF(ALAM.LT.XLAM(K)) THEN
HPLUS=HINP(K-1)+(HINP(K)-HINP(K-1))/(XLAM(K)-XLAM(K-1))*
* (ALAM-XLAM(K-1))
HOUT(I)=HOUT(I)+HPLUS*G(J)
ELSE IF(ALAM.EQ.XLAM(K)) THEN
HOUT(I)=HOUT(I)+HINP(K)*G(J)
ELSE
K=K+1
GO TO 130
END IF
200 CONTINUE
C
DO 220 J=1,NR1
A2=(J-1)*DLAM
ALAM=YLAM(I)-A2
K=K0
210 CONTINUE
IF(ALAM.GT.XLAM(K)) THEN
HMINUS=HINP(K)+(HINP(K+1)-HINP(K))/(XLAM(K+1)-XLAM(K))*
* (ALAM-XLAM(K))
HOUT(I)=HOUT(I)+HMINUS*G(J)
ELSE IF(ALAM.EQ.XLAM(K)) THEN
HOUT(I)=HOUT(I)+HINP(K)*G(J)
ELSE
K=K-1
GO TO 210
END IF
220 CONTINUE
IF(K0.LE.0) K0=1
300 CONTINUE
C
C 3. points near the end of the interval
C
IF(INTR1.LT.NLAMY) THEN
IF(MODE.EQ.1.AND.DLAM.GT.DLROT) THEN
SLAM=YLAM(NLAMY)
CALL KERNEL(MODE,NR,VROT,SLAM,FWHM,XLMAX,G)
DLAM=XLMAX/NR
END IF
K0=NLAMX
DO 500 I=IEND1,INTR1,-1
HOUT(I)=0.
DO 410 K=K0,1,-1
K0=K
IF(XLAM(K).LT.YLAM(I)) GO TO 420
410 CONTINUE
420 CONTINUE
DO 450 J=1,NR1
A2=(J-1)*DLAM
ALAM=YLAM(I)+A2
K=K0+1
430 CONTINUE
IF(ALAM.LT.XLAM(K)) THEN
HPLUS=HINP(K-1)+(HINP(K)-HINP(K-1))/(XLAM(K)-XLAM(K-1))*
* (ALAM-XLAM(K-1))
HOUT(I)=HOUT(I)+HPLUS*G(J)
ELSE IF(ALAM.EQ.XLAM(K)) THEN
HOUT(I)=HOUT(I)+HINP(K)*G(J)
ELSE
K=K+1
IF(K.GT.NLAMX) THEN
HOUT(I)=HOUT(I)+HP0*G(J)
GO TO 450
END IF
GO TO 430
END IF
450 CONTINUE
C
DO 470 J=1,NR1
A2=(J-1)*DLAM
ALAM=YLAM(I)-A2
K=K0
460 CONTINUE
IF(ALAM.GT.XLAM(K)) THEN
HMINUS=HINP(K)+(HINP(K+1)-HINP(K))/(XLAM(K+1)-XLAM(K))*
* (ALAM-XLAM(K))
HOUT(I)=HOUT(I)+HMINUS*G(J)
ELSE IF(ALAM.EQ.XLAM(K)) THEN
HOUT(I)=HOUT(I)+HINP(K)*G(J)
ELSE
K=K-1
GO TO 460
END IF
470 CONTINUE
IF(K0.LE.0) K0=1
500 CONTINUE
END IF
RETURN
END
C
C
C ****************************************************************
C
C
SUBROUTINE KERNEL(MODE,NR,VROT,SLAM,FWHM,XLMAX,G)
C -------------------------------------------------
C
C Kernel function for the rotational and/or instrumental convolution
C
C Input:
C
C MODE - for MODE = 1 - rotational convolution;
C for MODE = 2 - instrumental convolution,
C with a Gaussian instrumental profile
C NR - number of integration points for evaluating the
C convolution integrals
C VROT - v sin i [in km/s]
C SLAM - standard wavelength (for rotational convolution;
C no meaning for the instrumental convolution)
C FWHM - full width at half maximum of the instrum. profile function
C
C Output:
C
C XLMAX- width of the kernel
C G - normalized kernel function
C
implicit real*8 (a-h,o-z)
PARAMETER (MCONV=5000, MCONV1=MCONV+1)
PARAMETER (ONE=1.,
* TWO=2.,
* HALF=0.5)
PARAMETER (EPS=0.6,
* GAUSLM=3.,
* DLROT=10.)
DIMENSION G(MCONV1),gm(21)
c
data gm/ 1.128,
* 0.939, 0.773, 0.628, 0.504, 0.399,
* 0.312, 0.240, 0.182, 0.133, 0.101,
* 0.070, 0.052, 0.037, 0.024, 0.017,
* 0.012, 0.010, 0.009, 0.007, 0.006/
C
C set up integration limits for the convolution integral
C
IF(MODE.EQ.1) THEN
XLMAX=SLAM*VROT/2.997925E5
E1=TWO*(ONE-EPS)
E2=1.5707964*EPS
E3=ONE/3.1415927/XLMAX/(ONE-EPS/3.)
ELSE IF(MODE.EQ.2) THEN
XLMAX=GAUSLM*FWHM
D0=0.60056*FWHM
D1=0.5641896/D0
else if(mode.eq.3) then
vmac=vrot
xlmax=slam*vmac/2.997925e5*2.
nr=20
END IF
C
C evaluation of the kernel function G
C
DLAM=XLMAX/NR
NR1=NR+1
DO 10 J=1,NR1
X=(J-1)*DLAM
IF(MODE.EQ.1) THEN
X1=ABS(1.-(X/XLMAX)**2)
G(J)=(E1*SQRT(X1)+E2*X1)*E3
ELSE IF(MODE.EQ.2) THEN
G(J)=D1*EXP(-(X/D0)**2)
ELSE IF(MODE.EQ.3) THEN
G(J)=gm(j)
END IF
10 CONTINUE
C
C renormalization in order to have integral(G) = 1
C
SUM=0.
DO 20 J=2,NR
SUM=SUM+TWO*G(J)
20 CONTINUE
SUM=SUM+G(1)+G(NR1)
SUM=ONE/DLAM/SUM
C
DO 30 J=1,NR1
G(J)=G(J)*SUM
30 CONTINUE
c
c multiply by integration weights for trapezoidal integration
c
DO 40 J=1,NR1
G(J)=G(J)*DLAM
40 CONTINUE
G(1)=HALF*G(1)
G(NR1)=HALF*G(NR1)
RETURN
END
C
C
C ****************************************************************
C
C
SUBROUTINE INTERP(X,Y,XX,YY,NX,NXX,NPOL,ILOGX,ILOGY)
C ====================================================
C
C General interpolation procedure of the (NPOL-1)-th order
C
C for ILOGX = 1 logarithmic interpolation in X
C for ILOGY = 1 logarithmic interpolation in Y
C
C Input:
C X - array of original x-coordinates
C Y - array of corresponding functional values Y=y(X)
C NX - number of elements in arrays X or Y
C XX - array of new x-coordinates (to which is to be
C interpolated
C NXX - number of elements in array XX
C Output:
C YY - interpolated functional values YY=y(XX)
C
implicit real*8 (a-h,o-z)
DIMENSION X(1),Y(1),XX(1),YY(1)
EXP10(X0)=EXP(X0*2.30258509299405D0)
IF(NPOL.LE.0.OR.NX.LE.0) GO TO 200
IF(ILOGX.EQ.0) GO TO 30
DO 10 I=1,NX
10 X(I)=LOG10(X(I))
DO 20 I=1,NXX
20 XX(I)=LOG10(XX(I))
30 IF(ILOGY.EQ.0) GO TO 50
DO 40 I=1,NX
40 Y(I)=LOG10(Y(I))
50 NM=(NPOL+1)/2
NM1=NM+1
NUP=NX+NM1-NPOL
DO 100 ID=1,NXX
XXX=XX(ID)
DO 60 I=NM1,NUP
IF(XXX.LE.X(I)) GO TO 70
60 CONTINUE
I=NUP
70 J=I-NM
JJ=J+NPOL-1
YYY=0.
DO 90 K=J,JJ
T=1.
DO 80 M=J,JJ
IF(K.EQ.M) GO TO 80
T=T*(XXX-X(M))/(X(K)-X(M))
80 CONTINUE
90 YYY=Y(K)*T+YYY
YY(ID)=YYY
100 CONTINUE
IF(ILOGX.EQ.0) GO TO 130
DO 110 I=1,NX
110 X(I)=EXP10(X(I))
DO 120 I=1,NXX
120 XX(I)=EXP10(XX(I))
130 IF(ILOGY.EQ.0) RETURN
DO 140 I=1,NX
140 Y(I)=EXP10(Y(I))
DO 150 I=1,NXX
150 YY(I)=EXP10(YY(I))
RETURN
200 N=NX
IF(NXX.GE.NX) N=NXX
DO 210 I=1,N
XX(I)=X(I)
210 YY(I)=Y(I)
RETURN
END
C
C ********************************************************************
C ********************************************************************
synspec/synspec54.f
+23917
View File
File diff suppressed because it is too large Load Diff