;HW4: Fitting and Plotting a Spectrum
fits_read,'M7.fits',flux,header
crpos=where(strmid(header,0,6) eq "CRVAL1")
cr=strmid(header[crpos[0]],9,60)+0.
cdeltpos=where(strmid(header,0,6) eq "CDELT1")
cdelt=strmid(header[cdeltpos[0]],9,60)+0.
x=findgen(n_elements(flux))*cdelt+cr
!p.multi=[0,1,2]
c3=where(x ge 1750 AND x le 2050)
plot,x[c3],flux[c3],xtitle="Wavelength",ytitle="Flux",xrange=[1750,2050],xstyle=1
wcont=where((x ge 1800 AND x le 1820) OR (x ge 2020 AND x le 2040))
fit=linfit(x[wcont],flux[wcont])
oplot,x[c3],fit[0]+fit[1]*x[c3],linestyle=2
xyouts,1760,0.25,'Dashed line: continuum fit',charsize=1.5,/data
nor=flux/(fit[0]+fit[1]*x)
plot,x[c3],nor[c3],xtitle="Wavelength",ytitle="Normalized Flux",xrange=[1750,2050],xstyle=1
oplot,[1750,2050],[1,1]
xyouts,0.12,0.12,'Solid line: continuum',charsize=1.5,/normal
write_jpeg,'HW4.jpeg',tvrd()
header[crpos]='CRVAL1 =   1750.0'
fits_write,'HW4.fits',nor[c3],header
!p.multi=[0,1,1]
set_plot,'ps'
device,filename='HW4.ps'
plot,x[c3],nor[c3],xtitle="Wavelength",ytitle="Normalized Flux",xrange=[1750,2050],xstyle=1
oplot,[1750,2050],[1,1]
xyouts,1760,0.25,'Solid line: continuum',charsize=1.5,/data
device,/close
end