From e883e0cddf12d61d8fec3203788910ca6578f2fa Mon Sep 17 00:00:00 2001 From: Jason Eastman Date: Mon, 15 May 2017 15:57:29 -0400 Subject: [PATCH] adding more pro files --- exofast_plotchains.pro | 145 ++++++++++++++++++++++++++++ gettc.pro | 108 +++++++++++++++++++++ gridpars.pro | 60 ++++++++++++ mutual_eclipse.pro | 23 +++++ omc.pro | 184 ++++++++++++++++++++++++++++++++++++ recalc.pro | 23 +++++ resume_exofastv2.pro | 28 ++++++ sed/convertsed.pro | 65 +++++++++++++ sed/nextgenfin/getfiles.pro | 51 ++++++++++ 9 files changed, 687 insertions(+) create mode 100644 exofast_plotchains.pro create mode 100644 gettc.pro create mode 100644 gridpars.pro create mode 100644 mutual_eclipse.pro create mode 100644 omc.pro create mode 100644 recalc.pro create mode 100644 resume_exofastv2.pro create mode 100644 sed/convertsed.pro create mode 100644 sed/nextgenfin/getfiles.pro diff --git a/exofast_plotchains.pro b/exofast_plotchains.pro new file mode 100644 index 000000000..9573175a8 --- /dev/null +++ b/exofast_plotchains.pro @@ -0,0 +1,145 @@ +pro exofast_plotchains, ss, chainfile=chainfile + +if n_elements(chainfile) eq 0 then chainfile = 'chain.ps' + +burnndx = ss.burnndx +nsteps = ss.nsteps/ss.nchains +medndx = (nsteps-ss.burnndx)/2 +halfsigma = erf(1d/sqrt(2d0))/2d +lowsigndx = round(((nsteps-ss.burnndx)/2.d0 - (nsteps-ss.burnndx)*halfsigma)*ss.nchains) +hisigndx = round(((nsteps-ss.burnndx)/2.d0 + (nsteps-ss.burnndx)*halfsigma)*ss.nchains) +;; prepare the postscript device +mydevice=!d.name +set_plot, 'PS' +aspect_ratio=1.5 +xsize = 18 +ysize=xsize/aspect_ratio +!p.font=0 +device, filename=chainfile +device, set_font='Times',/tt_font +device, xsize=xsize,ysize=ysize,/color,bits=24 +charsize=1.5 +loadct,39,/silent +red = 254 + +allpars = [] +parnames = [] +medianpars = [] +!p.multi=[0,2,4] ;; 8 to a page + +chi2 = reform((*ss.chi2),ss.nsteps/ss.nchains,ss.nchains) +ymin = min(chi2[ss.burnndx:-1,*],max=ymax) +xtitle='!3 Chain link number' +if ymin lt 0 then sign = ' + ' else sign = ' - ' +ytitle='!3' + textoidl('\chi^2') + sign + strtrim(abs(ymin),2) +plot, [0],[0], psym=10, xtitle=xtitle, ytitle=ytitle,$ + charsize=charsize,xstyle=1, yrange=[0,ymax-ymin], xrange=[0,nsteps],font=1 +for l=0L, ss.nchains-1L do $ + oplot, chi2[*,l] - ymin, color=l*255d0/ss.nchains ;, transparency=100d0-100d0/ss.nchains +oplot, [ss.burnndx,ss.burnndx],[-9d99,9d99] +print, ymin + +for i=0, n_tags(ss)-1 do begin + for j=0, n_elements(ss.(i))-1 do begin + for k=0, n_tags(ss.(i)[j])-1 do begin + + derive = 0B + m=-1 + ;; this captures the detrending variables + if (size(ss.(i)[j].(k)))[1] eq 10 then begin + if (size(ss.(i)[j].(k)))[0] ne 0 then begin + for l=0L, n_tags(*(ss.(i)[j].(k)))-1 do begin + if (size((*(ss.(i)[j].(k))).(l)))[2] eq 8 then begin + for m=0L, n_elements((*(ss.(i)[j].(k))).(l))-1 do begin + if tag_exist((*(ss.(i)[j].(k))).(l)[m],'derive') then begin + if (*(ss.(i)[j].(k))).(l)[m].derive then begin + pars = (reform((*(ss.(i)[j].(k))).(l)[m].value,nsteps,ss.nchains)) +; pars = (*(ss.(i)[j].(k))).(l)[m].value[burnndx:*,*] + derive = 1B + label = (*(ss.(i)[j].(k))).(l)[m].label + unit = (*(ss.(i)[j].(k))).(l)[m].unit + latex = (*(ss.(i)[j].(k))).(l)[m].latex + best = (*(ss.(i)[j].(k))).(l)[m].best + endif + endif + endfor + endif + endfor + endif + endif else if n_tags(ss.(i)[j].(k)) ne 0 then begin + if n_tags(ss.(i)[j].(k)) ne 0 then begin + if tag_exist(ss.(i)[j].(k),'derive') then begin + if ss.(i)[j].(k).derive then begin + pars = (reform(ss.(i)[j].(k).value,nsteps,ss.nchains)) + derive = 1B + label = ss.(i)[j].(k).label + unit = ss.(i)[j].(k).unit + latex = ss.(i)[j].(k).latex + best = ss.(i)[j].(k).best + endif + endif + endif + endif + + if derive then begin + + if n_elements(allpars) eq 0 then allpars = transpose(pars) $ + else allpars = [allpars,transpose(pars)] + + ;; check for bad values + bad = where(~finite(pars),complement=good) + if bad[0] ne -1 then message, $ + "ERROR: NaNs in " + label + " distribution" + + ;; if angular, center distribution about the mode + if unit eq 'DEGREES' then halfrange = 180d0 $ + else if unit eq 'RADIANS' then halfrange = !dpi + if unit eq 'DEGREES' or unit eq 'RADIANS' then begin + + ;; find the mode + hist = histogram(pars,nbins=100,locations=x,/nan) + max = max(hist,modendx) + mode = x[modendx] + + toohigh = where(pars gt (mode + halfrange)) + if toohigh[0] ne -1 then pars[toohigh] -= 2.d0*halfrange + + toolow = where(pars lt (mode - halfrange)) + if toolow[0] ne -1 then pars[toolow] += 2.d0*halfrange + endif + + ;; 68% confidence interval + sorted = sort(pars) + medvalue = pars[sorted[medndx]] + upper = pars[sorted[hisigndx]] - medvalue + lower = medvalue - pars[sorted[lowsigndx]] + + if n_elements(medianpars) eq 0 then medianpars = [medvalue,upper,lower] $ + else medianpars = [[medianpars],[[medvalue,upper,lower]]] + + xmax = (medvalue + 4*upper) < max(pars) + xmin = (medvalue + 4*lower) > min(pars) + parnames = [parnames,latex] + + if xmin eq xmax then begin + message, 'WARNING: ' + label + ' is singularly valued.',/continue + endif else begin + + ;; plot labels + xtitle='!3 Chain link number' + ytitle='!3' + textoidl(latex + '_{' + ss.(i)[j].label + '}') + + plot, [0],[0], psym=10, xtitle=xtitle, ytitle=ytitle,$ + charsize=charsize,xstyle=1, yrange=[xmin,xmax], xrange=[0,nsteps],font=1 + + for l=0L, ss.nchains-1L do $ + oplot, pars[*,l], color=l*255d0/ss.nchains ;, transparency=100d0-100d0/ss.nchains + oplot, [ss.burnndx,ss.burnndx],[-9d99,9d99] + + endelse + endif + endfor + endfor +endfor + +end diff --git a/gettc.pro b/gettc.pro new file mode 100644 index 000000000..4a3efddbe --- /dev/null +++ b/gettc.pro @@ -0,0 +1,108 @@ +pro gettc, secondary=secondary + +npoints = 1d5 + +G = 2942.71377d0 ;; R_sun^3/(m_sun*day^2), Torres 2010 +nplanets = 3d0 +period = [1d0,2d0,10d0] +period = 1+randomu(seed,3)*30d0 +period = period[sort(period)] +tc = [0d0,3d0,5d0] +tc = randomu(seed,3)*max(period) +inc=[!dpi/2d0,!dpi/2d0,!dpi/2d0] +m1 = 1d0 +m2 = [0.001d0,0.01d0,0.01d0] +q = m1/m2 +a = (period^2*G*(m1+m2)/(4d0*!dpi^2))^(1d0/3d0) +;e = [0.5d0,0.25d0,0.9d0] +e = dblarr(3);[0.5d0,0.25d0,0.9d0] +omega = [!dpi/5d0,!dpi/6d0,!dpi/7d0] + + +bjd = min(tc) + dindgen(npoints)/(npoints-1)*max(period)*10 +tp = dblarr(nplanets) +for i=0, nplanets-1L do begin + phase = exofast_getphase(e[i],omega[i], /primary) + tp[i] = tc[i] - period[i]*phase +endfor + +if 0 then begin +b = exofast_getb2(bjd, inc=inc,a=a,tperiastron=tp, period=period,e=e,omega=omega,q=q,x0=x,y0=y,z0=z,x1=x1,y1=y1,z1=z1,x2=x2,y2=y2,z2=z2) + +for i=0, nplanets-1 do begin + + ndx = lclxtrem(b[i,*]) +!p.multi=[0,2,2] + plot, b[i,*] + oplot, b[i,ndx], b[i,ndx], psym=1, color='0000ff'x + +endfor +endif + +a2 = a*q/(1d0+q) +a1 = a2/q +tol = 1d-8 ;; 1 ms +ntransits = 10 + +for i=0L, nplanets-1 do begin + + print + print, 'planet ' + strtrim(i,2) + + for k=0, ntransits-1 do begin + + t0 = systime(/seconds) + minbjd = tc[i]-period[i]/4d0+k*period[i] + maxbjd = tc[i]+period[i]/4d0+k*period[i] + repeat begin + + bjd = (minbjd + maxbjd)/2d0 + + x1 = 0d0 + z1 = 0d0 + for j=0L, nplanets-1 do begin + meananom = (2.d0*!dpi*(1.d0 + (bjd - tp[j])/period[j])) mod (2.d0*!dpi) + eccanom = exofast_keplereq(meananom, e[j]) + trueanom = 2d0*atan(sqrt((1d0 + e[j])/(1d0 - e[j]))*tan(0.5d0*eccanom)) + + ;; calculate the star position in the barycentric frame + r1 = a1[j]*(1d0-e[j]^2)/(1d0+e[j]*cos(trueanom)) + + ;; rotate to observer's plane of reference + x1 += -r1*cos(trueanom + omega[j] + !dpi) + z1 += r1*sin(trueanom + omega[j] + !dpi)*sin(inc[j]) + endfor + + meananom = (2.d0*!dpi*(1.d0 + (bjd - tp[i])/period[i])) mod (2.d0*!dpi) + eccanom = exofast_keplereq(meananom, e[i]) + trueanom = 2d0*atan(sqrt((1d0 + e[i])/(1d0 - e[i]))*tan(0.5d0*eccanom)) + + r = a2[i]*(1d0-e[i]^2)/(1d0+e[i]*cos(trueanom)) + x = -r*cos(trueanom + omega[i])-x1 + z = r*sin(trueanom + omega[i] + !dpi)*sin(inc[i]) -z1 + + ;; bracket the primary or secondary + if keyword_set(secondary) then begin + if z le 0d0 and x gt 0d0 then minbjd = bjd $ + else if z le 0d0 and x le 0d0 then maxbjd = bjd $ + else if z gt 0d0 and x gt 0d0 then maxbjd = bjd $ + else if z gt 0d0 and x le 0d0 then minbjd = bjd + endif else begin + if z le 0d0 and x gt 0d0 then maxbjd = bjd $ + else if z le 0d0 and x le 0d0 then minbjd = bjd $ + else if z gt 0d0 and x gt 0d0 then minbjd = bjd $ + else if z gt 0d0 and x le 0d0 then maxbjd = bjd + endelse + + if x gt 0 then begin + maxbjd = bjd + endif else minbjd = bjd + + endrep until abs(x) le tol + + print, bjd, (bjd-tc[i]-k*period[i])*86400d0, systime(/seconds)-t0 + endfor +endfor + + +end diff --git a/gridpars.pro b/gridpars.pro new file mode 100644 index 000000000..03f9b16a6 --- /dev/null +++ b/gridpars.pro @@ -0,0 +1,60 @@ +pro gridpars + +files = file_search('sed/nextgenfin/nextgenfin2/*.idl',count=nfiles) + + +teff = [] +logg = [] +feh = [] +alpha = [] + +for i=0L, nfiles-1 do begin + + filename = file_basename(files[i],'.NextGen.spec.idl') + len = strlen(filename) + teff = [teff,double(strmid(filename,3,len-15))*100d0] + logg = [logg,double(strmid(filename,len-12,4))] + feh = [feh,double(strmid(filename,len-8,4))] + alpha = [alpha,double(strmid(filename,len-4,4))] + +endfor + +tefforig = teff +loggorig = logg +fehorig = feh +alphaorig = alpha + +teff = teff[sort(teff)] +teffs = teff[uniq(teff)] + +logg = logg[sort(logg)] +loggs = logg[uniq(logg)] + +feh = feh[sort(feh)] +fehs = feh[uniq(feh)] + +alpha = alpha[sort(alpha)] +alphas = alpha[uniq(alpha)] + +nteffs = n_elements(teffs) +nloggs = n_elements(loggs) +nfehs = n_elements(fehs) +nalphas = n_elements(alphas) +grid = bytarr(nteffs,nloggs,nfehs,nalphas) +for i=0,nteffs-1L do begin + for j=0, nloggs-1L do begin + for k=0, nfehs-1L do begin + for l=0, nalphas-1L do begin + match = where(tefforig eq teffs[i] and loggorig eq loggs[j] and $ + fehorig eq fehs[k] and alphaorig eq alphas[l]) + if match[0] ne -1 then begin + grid[i,j,k,l] = 1B + oplot, teffs[i],loggs[j] + endif + endfor + endfor + endfor +endfor + +stop +end diff --git a/mutual_eclipse.pro b/mutual_eclipse.pro new file mode 100644 index 000000000..4e74c5227 --- /dev/null +++ b/mutual_eclipse.pro @@ -0,0 +1,23 @@ +;; (x,y,z) = (0,0,0) is the center of the star +;; x -- An NPlanets x NTimes array of X positions of each planet, in stellar radii +;; y -- An NPlanets x NTimes array of Y positions of each planet, in stellar radii +;; z -- An NPlanets x NTimes array of Z positions of each planet, in stellar radii +;; p -- An NPlanets array of the radius of each planet, in stellar radii +;; u1 -- linear limb darkening parameter +;; u2 -- quadratic limb darkening parameter + +function mutual_eclipse, x, y, z, p, u1, u2 + +sz = size(x) +nplanets = sz[1] +ntimes = sz[2] +flux = dblarr(ntimes) + 1d0 + +for i=0, nplanets-1 do begin + + +endfor + +return, flux + +end diff --git a/omc.pro b/omc.pro new file mode 100644 index 000000000..861d5f76b --- /dev/null +++ b/omc.pro @@ -0,0 +1,184 @@ +;+ +; NAME: +; OMC +; +; PURPOSE: +; Calculates the best ephemeris from a list of transit times and +; plots the O-C diagram +; +; CALLING SEQUENCE: +; OMC, 'filename' [,period, t0=t0,/ps,epsname=epsname] +; +; INPUTS: +; +; FILENAME - The filename of a text file with columns T_C, error, +; and optionally telescope name. If telescope name is +; given, each unique telescope name will be plotted +; with a different color. The error must be in days. +; NOTE: "#" may be used as a comment string. +; +; OPTIONAL INPUTS: +; +; PERIOD - The period of the transit. If not given, the program +; will attempt to calculate it from the lowest common +; multiple of transits. +; NOTE: This may result in integer aliases depending on +; when the transits were taken. +; T0 - The zero point of the ephemeris. If not specified, the +; transit time closest to the error-weighted mean of all +; input times will be used, which minimizes the error in +; the ephemeris. +; EPSNAME - The name of the encapsulated postscript file; default +; is omc.eps. +; +; OPTIONAL KEYWORDS: +; +; PS - If set, an encapsulated postscript plot will be +; created. Otherwise, it will plot to the screen. It is +; not necessary to specify both PS and EPSNAME. +; +; OUTPUTS: +; +; 1) Plots the O-C diagram +; 2) Prints the best fit ephemeris with errors +; 3) Prints the chi^2/dof of the best fit +; 4) Prints the epoch, transit time, error (sec), residual (sec), +; and residual/err for each transit +; +; MODIFICATION HISTORY +; +; 2011/07 -- Written, Jason Eastman (OSU) +; +;- +pro omc, time, err, telescope=telescope, chi2=chi2, period=period, t0=t0, ps=ps, epsname=epsname + +;; read in the times +openr, lun, filename, /get_lun +line = "" +readf, lun, line +ncols = n_elements(strsplit(line)) +free_lun, lun +if ncols eq 2 then begin + readcol, filename, time, err, format='d,d', comment='#',/silent +endif else if ncols eq 3 then begin + readcol, filename, time,err,telescope, format='d,d,a', comment='#',/silent +endif else message, 'FILENAME must have 2 or three columns' + +;; Determine the lowest common multiple of transit times +;; this is either the period or an integer multiple of it +if n_elements(period) eq 0 then begin + niter = 1d0 + diff = (time-shift(time,1)) + good = where(diff gt 0.1) ;; if diff lt 0.1, it's the same transit + mindiff = min(diff[good]) + repeat begin + trialP = mindiff/niter + np = diff[good]/trialP + niter++ + if niter ge 10 then $ + message, 'Cannot determine best period; please input manually' + endrep until max(abs(np - round(np))) lt 0.1 + period = trialP +endif + +;; find the best t0 (i.e., nearest to the error-weighted mean of T_Cs) +if n_elements(t0) eq 0 then begin + tmean = total(time*err)/total(err) + np = round((time[0] - tmean)/period) + t0 = time[0] - np*period +endif + +;; calculate the epoch numbers (x axis) +epoch = round((time - t0)/period) + +;; fit a straight line to the Transit Times +;; TC = coeffs[0] +/- sigma[0] +;; P = coeffs[1] +/- sigma[1] +coeffs = poly_fit(epoch,time,1,measure_errors=err,yfit=yfit,sigma=sigma) + +;; calculate the residuals +omc = (time-yfit)*1440d0 +err *= 1440d0 ;; convert error to minutes + +;; prepare the plotting device (thanks astrobetter) +;; http://www.astrobetter.com/making-fonts-better-in-idl-postscript-output/ +if keyword_set(ps) or n_elements(epsname) ne 0 then begin + if n_elements(epsname) eq 0 then epsname = 'omc.eps' + mydevice = !d.name + set_plot, 'PS' + !p.font=0 + aspect_ratio=1.5 + xsize=9 + ysize=xsize/aspect_ratio + device, filename=epsname,/color,bits=24,encapsulated=1, /helvetica + device, xsize=xsize, ysize=ysize + LOADCT, 39,/silent + colors = [0,254,159,95,223,31,207,111,191,47] + charsizelegend = 0.09 + xlegend = 0.45 + ylegend = 0.50 + charsize = 0.5 +endif else begin + colors= ['ffffff'x,'0000ff'x,'00ff00'x,'ff0000'x,'0080ff'x,$ + '800080'x,'00ffff'x,'ffff00'x,'80d000'x,'660000'x] + charsizelegend = 0.03 + xlegend = 0.90 + ylegend = 0.99 + charsize = 1 +endelse +ncolors = n_elements(colors) + +;; plot the times, errors, and 0 line +ymin = min(omc-err) & ymax = max(omc+err) +xmin = min(epoch)-1 & xmax = max(epoch)+1 +xmin = -75 & xmax = 75 +plot, fix(epoch), omc,psym=3, yrange=[ymin,ymax],xrange=[xmin,xmax],$ + xtitle='Epoch',ytitle='O-C (min)', /xstyle;,xticks=4 +oploterr, epoch, omc, err, 3 +oplot, [-9d9,9d9],[0,0],linestyle=1 + +syms = [0,3,8,5,0,3,8,5] +fill = [1,1,1,1,0,0,0,0] +nsyms = n_elements(syms) + +;; overplot each telescope in different colors +if n_elements(telescope) ge 1 then begin + sorted = sort(telescope) + tnames = telescope[sorted[uniq(telescope[sorted])]] + for i=0, n_elements(tnames)-1 do begin + observed = where(telescopes eq tnames[i]) + if observed[0] ne -1 then begin + plotsym, syms[i mod nsyms], color=colors[i mod ncolors],fill=fill[i mod nsyms] + oplot, epoch[observed],omc[observed],psym=8 + xsize = (!x.crange[1] - !x.crange[0]) + ysize = (!y.crange[1] - !y.crange[0]) + xyouts, !x.crange[0] + xlegend*xsize,!y.crange[0]+(ylegend - i*charsizelegend)*ysize, $ + tnames[i],color=colors[i mod ncolors],charsize=charsize + oplot, [!x.crange[0]+xlegend*xsize-xsize/20],$ + [!y.crange[0]+(ylegend - (i-0.25)*charsizelegend)*ysize],psym=8 + endif + endfor +endif + +;; print the best-fit ephemeris, uncertainties, and chi^2 +print, 'T0 = '+string(coeffs[0],format='(f14.6)')+' +/- '+strtrim(sigma[0],2) +print, 'P = '+strtrim(coeffs[1],2) + ' +/- ' + strtrim(sigma[1],2) +print, 'chi^2 = '+strtrim(total(((omc)/err)^2),2) +print, 'dof = '+strtrim((n_elements(time)-2),2) +print, 'chi^2/dof = ' + strtrim(total(((omc)/err)^2)/(n_elements(time)-2),2) +print, '' + +;; print the epochs, times, errors, etc. +print, 'Epoch T_C ERR O-C (O-C)/err' +forprint, epoch, time, err*60, omc*60d0, omc/err,/t, $ + format='(i4,x,f14.6,x,i4,x,f7.2,x,f5.2)' + +;; close the postscript device +if keyword_set(ps) or n_elements(epsname) ne 0 then begin + device, /close + set_plot, mydevice + !p.font=1 +; spawn, 'gv ' + epsname + ' &' +endif + +end diff --git a/recalc.pro b/recalc.pro new file mode 100644 index 000000000..6b89e7a12 --- /dev/null +++ b/recalc.pro @@ -0,0 +1,23 @@ +pro recalc, idlfile + +if n_elements(idlfile) eq 0 then $;idlfile = 'EPIC229021605.c.mcmc.idl' + idlfile = 'GJ1132.c.noteffprior.demc2.mcmc.idl' +restore, idlfile +prefix='EPIC229021605.c2.' +;; output filenames +basename = file_basename(prefix) +label = "tab:" + basename +caption = "Median values and 68\% confidence interval for " + basename +parfile = prefix + 'pdf.ps' +covarfile = prefix + 'covar.ps' +texfile = prefix + 'median.tex' +chainfile = prefix + 'chain.ps' + +if keyword_set(covar) then nocovar = 0 $ +else nocovar = 1 + +exofast_plotdist2, mcmcss, pdfname=parfile, covarname=covarfile,nocovar=nocovar +exofast_latextab2, mcmcss, caption=caption, label=label,texfile=texfile +exofast_plotchains, mcmcss, chainfile=chainfile + +end diff --git a/resume_exofastv2.pro b/resume_exofastv2.pro new file mode 100644 index 000000000..35a8d51a0 --- /dev/null +++ b/resume_exofastv2.pro @@ -0,0 +1,28 @@ +;+ +; NAME: +; resume_exofastv2 +; PURPOSE: +; Given an IDL save file generated by exofastv2, resumes and extends +; the MCMC portion of the fit. +; +; DESCRIPTION: +; +; CALLING SEQUENCE: +; resume_exofastv2, 'savefile', maxsteps, nthin +; +; INPUTS: +; SAVEFILE - The filename of an IDL save file generated by +; EXOFASTv2. +; MAXSTEPS - The maximum number of steps to take. This is the +; total number of steps, not the number of additional steps. +; +; +; OPTIONAL INPUTS: +; +pro resume_exofastv2, savfile, maxsteps + +restore, savfile + + + +end diff --git a/sed/convertsed.pro b/sed/convertsed.pro new file mode 100644 index 000000000..e0da2a63b --- /dev/null +++ b/sed/convertsed.pro @@ -0,0 +1,65 @@ +;; this takes the raw model flux files (retrieved with getfiles.pro) +;; into the required format for exofast_interp_model.pro. Restricting +;; these files to 0.1 to 24 microns reduces the file sizes, load +;; times, and RAM footprint dramatically. Interpolating now saves us +;; time later but may lose generality +pro convertsed + +files = file_search(['/home/jeastman/nextgenfin/nextgen.?????.spec',$ + '/h/onion0/idl/EXOFASTv2/sed/nextgenfin/nextgenfin2/nextgen.?????.spec'],$ + count=nfiles) + +w1 = findgen(24000)/1000+0.1 ;; wavelength scale on which to interpolate (microns) + +line = '' +for i=0L, nfiles-1 do begin + + openr, lun, files[i], /get_lun + readf, lun, line ;; discard uselss header + readf, lun, line + teff = double((strsplit(line,/extract))[3]) + readf, lun, line + logg = double((strsplit(line,/extract))[3]) + readf, lun, line + meta = double((strsplit(line,/extract))[3]) + readf, lun, line + alpha = double((strsplit(line,/extract))[3]) + free_lun,lun + + teffstr = strtrim(round(teff/100.0),2) + + if logg lt 0 then loggsign = '-'$ + else loggsign = '+' + + if meta lt 0 then zsign = '-'$ + else zsign = '+' + + if alpha lt 0 then asign = '-'$ + else asign = '+' + + fmt = '("lte",a,a,f3.1,a,f3.1,a,f3.1,".NextGen.spec.idl")' + outname = string(teffstr, loggsign, abs(logg), zsign, abs(meta), asign, abs(alpha), format=fmt) + + if not file_test(outname) then begin + readcol, files[i], wavelength, flux, comment='#', format='f,f',/silent + + lamflam=flux*wavelength + lamflam1=interpol(lamflam,wavelength,w1*1d4) ;; erg/s/cm^2 + + save, lamflam1, filename=outname + +; openw, lun, outname, /get_lun +; printf, lun, string(teff,logg,meta,alpha,format='(i,x,f4.1,x,f4.1,x,f4.1," Teff, logg, [M/H], alpha")') +; printf, lun, strtrim(24000,2) + ' number of wavelength points. lambda = findgen(24000)/1000+0.1 ;; microns' +; printf, lun, lamflam1 +; free_lun, lun + + endif + + print, 'done with ' + outname + ' (' + files[i] + ')' + +endfor + +stop + +end diff --git a/sed/nextgenfin/getfiles.pro b/sed/nextgenfin/getfiles.pro new file mode 100644 index 000000000..10dcf8548 --- /dev/null +++ b/sed/nextgenfin/getfiles.pro @@ -0,0 +1,51 @@ +pro getfiles + + + +for i=1L, 13631L do begin + + spawn, 'curl "http://svo2.cab.inta-csic.es/theory/newov/ssap.php?model=bt-nextgen-agss2009&fid=' + strtrim(i,2) + '&format=ascii"', output + teff = (strsplit(output[1],/extract))[3] + logg = (strsplit(output[2],/extract))[3] + meta = (strsplit(output[3],/extract))[3] + alpha = (strsplit(output[4],/extract))[3] + + npoints = n_elements(output)-10 + wavelength = dblarr(npoints) + flux = dblarr(npoints) + + for j=0L, npoints-1 do begin + arr = strsplit(output[j+9],/extract) + wavelength[j] = arr[0] + flux[j] = arr[1] + endfor + + outname = string(round(double(teff)/100.0), double(logg), double(meta), format='("lte",i2,"-",f3.1,"-",f3.1,".NextGen.spec")') + + openw, lun, outname, /get_lun + printf, lun, teff, logg, meta, alpha 'Teff, logg, [M/H], alpha' + printf, lun, npoints, 'number of wavelength points' + printf, lun, wavelength + printf, flux + free_lun, lun + + stop + +endfor + + +;readcol, 'filelist.txt', files, format='a' +;for i=0L, nfiles-1 do begin +; +; extractedfile = (strsplit(files[i],'.gz',/extract,/regex))[0] +; if ~file_test(files[i]) and ~file_test(extractedfile) then begin +; arr = strsplit(files[i],'-.',/extract) +; zfolder = 'Z-' + arr[3] + '.' + arr[4] +; +; spawn, 'wget ftp://calvin.physast.uga.edu/pub/NextGen/Spectra/' + zfolder + '/' files[i] +; +; endif +; +;endfor + +end