pro sp_generate_model_OH_line_spectrum ; ; generates a Mimir-like spectrum from the ; known list of OH night sky lines ; ; DPC 20060304 original code ; filein = dialog_pickfile(Title="Select OH line list file") openr, lu_lines, filein, /GET_LUN ; ; device, decomposed=0 black = 14 loadct,0 tvlct, [0, 255, 0, 0, 255], [0, 0, 255, 0, 255], [0, 0, 0, 255, 255], 10 ; ; what spectral limits? ; print, "Enter starting, ending wavelengths [in microns] to model" read, wave_start, wave_stop print, "Enter wavelength FWHM resolution [in nm]" read, resol print, "Enter spectral channel width [in nm]" read, width ; ; read the OH line list ; OH_wave = fltarr(10000) OH_strength = fltarr(10000) iline = 0 ; while ~EOF(lu_lines) do begin line = '' readf, lu_lines, line test = STRPOS(line, "#") if(test ne 0) then begin wave = 0.0 strength = 0.0 reads, line, wave, strength ; ; keep only the stronger ones ; ;if(strength gt 10.) then begin OH_wave[iline] = wave/1.0e4 OH_strength[iline] = strength ;print, wave, strength iline++ ;endif endif endwhile FREE_LUN, lu_lines nlines = iline OH_wave = OH_wave[0:nlines-1] OH_strength = OH_strength[0:nlines-1] print, "read in ",nlines," OH line parameters" ; ;----------- lines all read in ---------------- ; ; next, convolve with gaussian of indicated resolution, ; for 10x number of channels desired ; nch = ROUND(((wave_stop - wave_start) * 1.0e4)/ width) + 101 print, "Nch = ",nch ; xwave = (indgen(nch) * width/1.0e4) + (wave_start - 0.005 * width) ; offset by 5 widths spec = dblarr(nch) ; ; for each wavelength, compute combined strength of all ; gaussians (one per OH line) ; which_lines = intarr(nlines) ; for ich = 0, nch-1 do begin spec[ich] = 0.0d0 for iline = 0, nlines-1 do begin ; ; if more than 100 sigma away, don't compute gaussian ; dwave = (OH_wave[iline] - xwave[ich]) test = double(abs(dwave/(resol*1.0e-3/2.354))) if(test lt 10.) then begin ; ; ok, within +/- 100 sigma, compute strength at this ; wavelength ; gauss = OH_strength[iline] * exp(-1.0d0*(test^2)/2.0d0) ;print, "Ich = ",ich," iline = ",iline," dwave = ",dwave," test = ",test," gauss = ",gauss spec[ich] = spec[ich] + gauss which_lines[iline] = iline endif endfor endfor indw = WHERE(which_lines ne 0,nwline) which_lines = which_lines[indw] print, "Found, modeled a total of ",nwline," OH lines in the desired range" window, 11, xsize=900, ysize=500, xpos=200, ypos=200, title="Smoothed High-Resolution Spectrum" ymin = min(spec) ymax = max(spec) xmin = min(xwave) xmax = max(xwave) ymplot = (ymax-ymin)*0.1+ymax plot, xwave, spec,xtitle="Wavelength [microns]", $ ytitle="OH line strength", xrange=[xmin,xmax],xstyle=1, color=black, $ yrange=[ymin, ymplot] ;cursor, xx, yy, /wait ;wait, 0.5 ; ; add vertical lines for all OH lines in range ; for iwline = 0, nwline-1 do begin iline = which_lines[iwline] xpline = OH_wave[iline] if(xpline ge xmin and xpline le xmax) then begin plots, [xpline,xpline],[0.,2.0*OH_strength[iline]], color=12 endif endfor ; ; add text about actions ; xmessage = 0.1*(xmax - xmin)+xmin ymessage = 0.97*(ymplot-ymin) + ymin str = STRING("Select OH line, or click outside box to exit") xyouts, xmessage, ymessage, str, color=11 ; another: wshow, 11 cursor, xx, yy, /wait, /data wait, 0.5 ; ; test for cursor in box ; if(xx ge xmin and xx le xmax and yy ge ymin and yy le ymplot) then begin ; ; cursor in box, find nearest peak in convolved spectrum ; and fit a gaussian to its center ; wave_dif = abs(xwave - xx) min_dif = MIN(wave_dif) ; if(N_ELEMENTS(min_dif) ne 1) then begin min_dif = min_dif[0] endif g_ind = WHERE(wave_dif eq min_dif) ; channel number of cursor ; ; grab nearest +/- 15 spectral channels (+/- FWHM) for plot ; chmin = 0 > (g_ind-15) chmax = (g_ind+15) < (nch-1) xgaus = xwave[chmin:chmax] ygaus = spec[chmin:chmax] est = [max(ygaus),xwave[g_ind],width/(2.354e3)] yfit = GAUSSFIT(xgaus, ygaus, a, NTERMS=3, SIGMA=s_a, ESTIMATES=est) ; ; plot yfit in red if strength > 0 ; if(a[0] gt 0.) then begin oplot, xgaus, yfit, color=11 ; ; and print center, uncertainty for center, peak ; print, format="('OH convolved line center = ',f7.5,' +/- ',f7.5,' microns. Strength = ',f9.1)",$ a[1],s_a[1],a[0] endif else print," Bad fit, try again" ; ; ; ;wave_dif = abs(OH_wave - xx) ;min_dif = MIN(wave_dif) ;if(N_ELEMENTS(min_dif) ne 1) then begin ; min_dif = min_dif[0] ;endif ;oh_ind = WHERE(wave_dif eq min_dif) ;x_OH = OH_wave[oh_ind] ;print,format="('OH line center = ',f8.6,' microns. Strength = ',f10.1)",x_OH, OH_strength[oh_ind] ; ; mark this line red ; ;plots, [x_OH,x_OH],[0.,2.0*OH_strength[oh_ind]], color=11 goto, another endif ; ; now, resample spectrum with desired channel width ; mch = ROUND(((wave_stop - wave_start) * 1.0e3)/ width) + 1 zwave = indgen(mch) * width / 1.0e3 + wave_start final_spec = dblarr(mch) print, "mch = ",mch ; ; perform block averaging ; for ich = 0, mch-1 do begin nch0 = (10*ich + 45) > 0 nch1 = nch0 + 9 for nch = nch0, nch1 do begin final_spec[ich] = final_spec[ich] + spec[nch] endfor endfor ; zmin = min(zwave) zmax = max(zwave) ymax = max(final_spec) ymin = min(final_spec) ypmax = (ymax-ymin)*0.1+ymax window, 10, xsize = 900, ysize=500, title="Spectrum Sampled at Indicated Channel Width", xpos=200, ypos=200 plot, zwave, final_spec, xtitle="Wavelength [microns]",ytitle="OH Sky Lines Strength", $ xrange=[zmin, zmax],yrange=[ymin, ypmax],xstyle=1, psym=10, color=black ; ; xmessage = 0.1*(zmax-zmin)+zmin ymessage = 0.97*(ypmax-ymin)+ymin str = "Click inside box to overlay unsmoothed spectrum, outside to exit" xyouts, xmessage, ymessage, str, color=11 ; overlay: wshow, 10 cursor, xx, yy, /wait, /data wait, 0.5 ; ; test for cursor in plot box ; if(xx ge zmin and xx le zmax and yy ge ymin and yy le ypmax) then begin ; ; overlay desired ; scale to same max ; ysmax = max(spec) ratio = ymax/ysmax oplot, xwave, spec*ratio, color=13 goto, overlay endif ; fileout = dialog_pickfile(title="Name of Output Spectrum File") openw, lu_out, fileout, /GET_LUN for ich = 0, mch-1 do begin printf, lu_out, format="(f12.8,1x,e15.6)",zwave[ich],final_spec[ich] endfor FREE_LUN, lu_out wdelete, 11 wdelete, 10 end