program lineextractor
implicit none

! ------------ DESCRIPTION ---------
!
! Extract line luminosities from a nebular spectrum
! as described in Jerkstrand et al 2012
!
! --------------------------------------------------

INTEGER :: i, ind, bin, Nspectra, j, m
DOUBLE PRECISION, DIMENSION(100000) :: wl, flux, fluxsel, time
INTEGER, PARAMETER :: Nlines=16
DOUBLE PRECISION, DIMENSION(Nlines) :: linelow, linehigh, linecenter, vcut, vcont, A_cont, B_cont, contr, contl, contwlr, contwll, linetot, moment1
DOUBLE PRECISION :: corrf, c, distance_mpc, vcut0, vcont0, df, scaling, use5577
CHARACTER :: trash
CHARACTER(len=200), DIMENSION(100) :: filein
DOUBLE PRECISION :: totlim, totlim_sel
DOUBLE PRECISION, DIMENSION(300) :: trashlong

totlim = 0
totlim_sel = 0

! INPUT ----------------------------------------

distance_mpc = 5.5  ! distance to SN in MPc

Nspectra = 1  ! Number of spectra to process

time(1:Nspectra) = (/100/) ! Epoch of spectra 

filein(1) = '/Users/anders/data/2012aw/332d/SN2012aw_56334.dat_photcal0.86_EBV0.074_V778'  ! Filenames

vcut0 = 2700  !  2700 km/s in Jerkstrand et al 2012 (Type IIP) 3500 km/s in Jerkstrand et al 2015 (Type IIb)

vcont0 = 1.25*vcut0  ! 1.25 is standard


! -END INPUT ------------------------------------

linetot = 0
moment1 = 0
contr = 100
contl = 100
wl = 0
flux = 0
ind = 0
linecenter = 0.

vcut(:) = vcut0
vcont(:) = vcont0

linelow(1) = 4571.
linehigh(1) = 4571.
linecenter(1) = 4571.

use5577 = 5577
linelow(2) = use5577
linehigh(2) = use5577
linecenter(2) = use5577

linelow(3) = 5890.
linehigh(3) = 5896
linecenter(3) = 5893.

linelow(4) = 6300.
linehigh(4) = 6364.
linecenter(4) = 6316.

linelow(5) = 6563.
linehigh(5) = 6563.
linecenter(5) = 6563.

linelow(6) = 7155.
linehigh(6) = 7155.
linecenter(6) = 7155.

linelow(7) = 7291.
linehigh(7) = 7323.
linecenter(7) = 7303.

linelow(8) = 7771.
linehigh(8) = 7775.
linecenter(8) = 7774.

linelow(9) = 8498.
linehigh(9) = 8662.

linelow(10) = 9260.  
linehigh(10) = 9266.  
linecenter(10) = 9263.

linelow(11) = 16000.
linehigh(11) = 16400.

linelow(12) = 15040
linehigh(12) = 15040

linelow(13) = 11286
linehigh(13) = 11302

linelow(14) = 13163
linehigh(14) = 13165

linelow(15) = 106822
linehigh(15) = 106822

linelow(16) = 66360
linehigh(16) = 66360


! ---------------------------------------------------------

c = 3e10

open(20,file='fluxes.dat',STATUS='REPLACE')

DO j = 1, Nspectra

   open(10,file=filein(j),STATUS='OLD')

   ind = 0; wl = 0; flux = 0; linetot = 0
   contr = 100; contl = 100; 

   DO 
      
      read(10,*,end=300), wl(ind+1), flux(ind+1) ! Read wavelength (A) and F_lambda (erg s^-1 cm^-2 A^-1)

      IF (flux(ind+1) < 0) CYCLE  ! go to next data point
 
      ind = ind + 1

      DO i = 1, Nlines

         IF (wl(ind) > linehigh(i) .AND. wl(ind) < linehigh(i)*(1 + vcont(i)/3e5) .AND. flux(ind) < contr(i))  THEN  ! Wavelength must be right of line center, within the continuum wavelength limit, and flux has to be lowest one so far

            contr(i) = flux(ind)

            contwlr(i) = wl(ind)

         END IF

         IF ( (wl(ind) < linelow(i)*(1-vcont(i)/3.0/3e5)) .AND. (wl(ind) > linelow(i)*(1 - vcont(i)/3e5)) .AND. (flux(ind) < contl(i)) .AND. (flux(ind) > 0)) THEN  !TEMP wavelength must be smaller than 1/3 of the way to the continuum limit : this is due to some problems in O I 5577, a crazy continuum is selected unless guided 
          
            contl(i) = flux(ind)

            contwll(i) = wl(ind)
        
         END IF

      END DO

   END DO

   300 close(10)

   DO i = 1, Nlines  ! Find continuum line, y = A + Bx

      B_cont(i) = (contr(i) - contl(i))/(contwlr(i) - contwll(i))

      A_cont(i) = contr(i) - B_cont(i)*contwlr(i)

      print *, 'Continuum', i, B_cont(i), A_cont(i), contwll(i), contwlr(i)

   END DO

   DO i = 2, ind

      DO m = 1, Nlines

         !IF ((wl(i) > linelow(m)*(1 - vcut(m)/3e5)) .AND. (wl(i) < linehigh(m)*(1. + vcut(m)/3e5)) .AND. (flux(i) > A_cont(m) + B_cont(m)*wl(i)) .and. (wl(i) > contwll(m)) .and. (wl(i) < contwlr(m))) THEN ! TEMP added wavelength constraints : DONT allow contributions from wavelengths outside where continuum has been determined

         IF ((wl(i) > linelow(m)*(1 - vcut(m)/3e5)) .AND. (wl(i) < linehigh(m)*(1. + vcut(m)/3e5))) THEN

            linetot(m) = linetot(m) + (flux(i)-(A_cont(m) + B_cont(m)*wl(i)))*(wl(i)-wl(i-1))

            moment1(m) = moment1(m) + (flux(i)-(A_cont(m) + B_cont(m)*wl(i)))*wl(i)*(wl(i)-wl(i-1))

         END IF

      END DO
       
   END DO
  
   df = 4*3.1415*(distance_mpc*1e6*3.086d18)**2

   DO i = 1, Nlines

      write(*,'(F10.1, 3ES12.2, F10.2, ES12.2, F10.2, 5ES12.4)') linelow(i), linetot(i)*df, A_cont(i), B_cont(i), contwll(i), contl(i), contwlr(i), contr(i), moment1(i)/linetot(i), (moment1(i)/linetot(i)-linecenter(i))/linecenter(i)*3e5  

   end do


END DO
 
END program fluxc
