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