program oh_lines
c+
c do various things with NIR OH spectrum file provided by
c Francois Piche
c
c
c input units are:
c angstroms and photons/s/m**2/arcsec**2
c
c 10751.0 49
c 10772.1 16
c 10831.3 60
c
c e=hv
c
c h=6.626 10^-34 J s
c =6.626 10^-27 erg s
c
c
c
c* The monochromatic flux reaching the Earth's atmosphere from Vega at
* 5556ang is
*
* 3.39**-9 ergs cm-2 s-1 A-1
* 3.50**-20 ergs cm-2 s-1 Hz-1
* 948 photons cm-2 s-1 A-1
*
* m(lam)-m(5556)=-2.5log(Flam(lam)/Flam(5556))
* therefore Flam(lam)= Flam(5556)*10**( /-2.5)
*
* Some useful cgs to SI conversions
*
* 1 Jy= 10**-23 ergs cm-2 Hz-1 = 10**-26 W m-2 Hz-1
* 1 W m-2 Hz-1 = 10**3 ergs sec-1 cm-2 Hz-1
*
* AB magnitudes defined by Oke and Gunn(1983)
* AB=-2.5log(Fnu) -48.60
* where 48.60 comes from 2.5log(Fnu(5556)) from above we get 48.63
* the difference is due to the different assumed zero-point for
* alpha-lyrae.
* The AB system is not actually zero-pointed on alpha-lyrae
* For convenience to that the numbers are similar it is
* zero-pointed on an assumed flux for alpha-lyrae which has since
* been superseded, if you know what I mean!
*
c-
implicit none
integer ioerr,irec
integer ibin,ibin_max
integer iline,nlines
integer iwave,nwave
integer iunit_out1,iunit_out2,iunit_out3
integer iunit_out4,iunit_out5,iunit_out6
real wave0,wave_final,delta_wave
real filter_width
real delta_wave_filter
real wave_microns(10000)
real wave_ang(10000)
real flux(10000)
real flux_cgs(10000)
real flux_si(10000)
real wave_rebinned(10000)
real flux_rebinned(10000)
real flux_rebinned_broad(10000)
real c
real flux_j,flux_h
real flux_j_cgs,flux_h_cgs
real flux_j_si,flux_h_si
real flux_ratio
character record*132
character filename_in*132
character filename_out*132
c=3e8 ! m s^-1
iunit_out1=6
iunit_out2=7
iunit_out3=8
iunit_out4=9
iunit_out5=10
iunit_out6=11
filename_in='oh_lines.txt'
open(1,file=filename_in)
write(*,*) filename_in(1:40),' opened '
ioerr=0
irec=0
write(*,*) ' irec,wave_ang(irec),flux(irec),flux_cgs(irec),flux_si(irec) '
dowhile(ioerr .eq. 0)
read(1,'(a)',iostat=ioerr) record
* write(*,*) ioerr,record
if(record(1:1) .ne. '#' .and. ioerr .eq. 0) then
irec=irec+1
read(record,*) wave_ang(irec),flux(irec)
* convert wavelengths to microns
wave_microns(irec)=wave_ang(irec)/10000.0
* convert to ergs/s/cm^2/arcsec**2
* input units are photons/s/m**2/arcsec**2
* but the wavelengths are stored internally as microns
*
* e=h nu
*
* h=6.626 10^-34 J s
* =6.626 10^-27 erg s
*
* convert to erg/s/cm^2/arcsec^2
flux_cgs(irec)=
: flux(irec)*(3e18/wave_ang(irec))*6.63e-27*1e-4
write(*,*) irec,wave_ang(irec),flux(irec),flux_cgs(irec),flux_si(irec)
endif
enddo
nlines=irec
write(*,*) nlines,' lines read in from oh_lines.txt '
* pass through the data and bin into low resolution ie 0.001microns
delta_wave=0.001 ! microns == 10ang
wave0=1.0
wave_final=2.50
nwave=(wave_final-wave0)/delta_wave
ibin_max=((wave_final-wave0)/delta_wave)+1
do iline=1,nwave
ibin=((wave_microns(iline)-wave0)/delta_wave)+1
flux_rebinned(ibin)=flux_rebinned(ibin)+flux_cgs(iline)
ibin_max=max(ibin_max,ibin)
enddo
write(*,*) ' ibin_max = ',ibin_max
* cycle through the bins and compute the wavelengths
do ibin=1,ibin_max
wave_rebinned(ibin)=wave0+((ibin-1)*delta_wave)
* write(*,*) ibin,wave_rebinned(ibin),flux_rebinned(ibin)
write(iunit_out3,*) wave_rebinned(ibin),flux_rebinned(ibin)
* write(2,*) ibin,wave_rebinned(ibin)-(delta_wave/2),' 0.0 '
* write(2,*) ibin,wave_rebinned(ibin)-(delta_wave/2),
* : flux_rebinned(ibin)
* write(2,*) ibin,wave_rebinned(ibin)+(delta_wave/2),
* : flux_rebinned(ibin)
* write(2,*) ibin,wave_rebinned(ibin)+(delta_wave/2),' 0.0 '
enddo
* compute the effective sky within a narrow band filter
* by adding up all the lines that lie within a filter
call getr4('Enter the percentage filter width : ',filter_width)
filter_width=filter_width/100.0
wave0=1.0
delta_wave=0.001
wave_final=2.50
nwave=(wave_final-wave0)/delta_wave
nwave=nwave+1
do iwave=1,nwave
wave_rebinned(iwave)=wave0+((iwave-1)*delta_wave)
flux_rebinned(iwave)=0
delta_wave_filter=wave_rebinned(iwave)*filter_width
* sum lines that would lie in filter
do iline=1,nlines
if(abs(wave_microns(iline)-wave_rebinned(iwave)) .lt.
: (delta_wave_filter/2.0)) then
flux_rebinned(iwave)=flux_rebinned(iwave)+flux_cgs(iline)
endif
enddo
write(iunit_out4,*) iwave,wave_rebinned(iwave),flux_rebinned(iwave)
enddo
write(*,*) nlines,' lines used '
write(*,*) ' wavelength range : ',wave_rebinned(1),wave_rebinned(nwave)
* compute the flux in J filter
flux_j_cgs=0
flux_j_si=0
do iline=1,nlines
if(wave_microns(iline) .ge. 1.00 .and.
: wave_microns(iline) .le. 1.40 ) then
flux_j_cgs=flux_j_cgs+flux_cgs(iline)
flux_j_si=flux_j_si+flux_si(iline)
write(*,*) iline,wave_microns(iline),flux_cgs(iline),flux_j_cgs
endif
enddo
* compute the flux in H filter
flux_h_cgs=0
flux_h_si=0
do iline=1,nlines
if(wave_microns(iline) .ge. 1.40 .and.
: wave_microns(iline) .le. 1.90 ) then
flux_h_cgs=flux_h_cgs+flux_cgs(iline)
flux_h_si=flux_h_si+flux_si(iline)
write(*,*) iline,wave_microns(iline),flux_cgs(iline),flux_h_cgs
endif
enddo
write(*,*) ' Integrated line flux in J : ',flux_j_cgs,flux_j_si
write(*,*) ' Integrated line flux in H : ',flux_h_cgs,flux_h_si
* compute the relative flux in the narrow band filter relative to a broad
* band filter
* J
do iwave=1,nwave
if(wave_rebinned(iwave) .ge. 1.10 .and.
: wave_rebinned(iwave) .le. 1.35) then
flux_ratio=1e4
if(flux_rebinned(iwave) .gt. 0) flux_ratio=flux_j_cgs/flux_rebinned(iwave)
write(iunit_out5,*) wave_rebinned(iwave),flux_rebinned(iwave),
: flux_ratio
write(6,*) wave_rebinned(iwave),flux_rebinned(iwave),flux_j_cgs,
: flux_ratio
endif
enddo
* H
do iwave=1,nwave
if(wave_rebinned(iwave) .ge. 1.50 .and.
: wave_rebinned(iwave) .le. 1.80) then
flux_ratio=1e4
if(flux_rebinned(iwave) .gt. 0) flux_ratio=flux_j_cgs/flux_rebinned(iwave)
write(iunit_out5,*) wave_rebinned(iwave),flux_rebinned(iwave),
: flux_ratio
write(6,*) wave_rebinned(iwave),flux_rebinned(iwave),flux_j_cgs,
: flux_ratio
endif
enddo
end