c     ***********************************************************
c     *                                                         *
c     *           CELIAS/SEM analysis software                  *
c     *                 March 2000                              *
c     *           USC Space Sciences Center                     *
c     *                                                         *
c     * This program converts CELIAS/SEM 5 minute averaged data *
c     * into 26-34 nm and 0.1-50 nm flux values.                *
c     *                                                         * 
c     * The input files are (1) qefficiency.dat, (2) flux.dat,  *
c     * (3) trans.dat, (4) background.dat and (5) carbon.dat.   *
c     *                                                         *
c     * The output file is named semflux.dat.                   *
c     *                                                         *
c     *          Written by: Pradip Gangopadhyay                *
c     *          USC Space Sciences Center                      *
c     *                                                         *
c     *          Contact: spacectr@usc.edu                      *
c     *                                                         *
c     ***********************************************************
      dimension eff(50),flux(50),trans(50),coeff(50)
      dimension ch1a(20000),ch1b(20000),ch2(20000)
      dimension time(20000),fl1(20000),fluxn(50),fl1old(20000)
      dimension carb(20000),q1(15),q3(15),sol(15),xmu(15)
      dimension beta1(15),beta3(15),r(20000)
      dimension fl3old(20000),fl3new(20000),fl1new(20000)
      dimension sd1(20000),sd2(20000),sd3(20000),x1(20000)
      dimension y1(20000),z1(20000),r1(20000)
      integer*4 ii,isave,j,lk,iyr,iday(20000),jsec(20000)
      character filename*50
c
c     ***********************************************************
c     itot is the number of 10 A bins between 0 and 500 Angstroms
c     area is the SEM aperture area in cm**2
c     ***********************************************************
c
      itot=50
      area=0.204362
c
c     ***********************************************************
c     store the quantum efficiency of the three channels q1(i), 
c     eff(i) and q3(i) for channels 1, 2 and 3 respectively.
c     ***********************************************************
c
      open(1,file='qefficiency.dat',status='old')
      do 10 i=1,itot
         read(1,*)ia,eff1,eff(i),eff3
         if(i.eq.18)q1(i-17)=eff1
         if(i.eq.18)q3(i-17)=eff3
         if((i.gt.25).and.(i.lt.36))then
            q1(i-24)=eff1
            q3(i-24)=eff3
         end if
 10   continue
      close(1)
c
c     ************************************************************
c     store the composite solar flux spectrum needed for flux.
c     flux(i) is used for the central channel flux while sol(i) is
c     the solar flux used for channels 1 and 3.
c     ************************************************************
c
      open(2,file='flux.dat',status='old')
      do 20 i=1,itot      
         read(2,*)a,flux(i)
         if(i.eq.18)sol(i-17)=flux(i)
         if((i.gt.25).and.(i.lt.36))sol(i-24)=flux(i)
 20   continue
      close(2)
c
c     ************************************************************
c     store the coefficient of the carbon transmission coeff(i).
c     ************************************************************
c
      open(3,file='trans.dat',status='old')
      do 30 i=1,itot
         read(3,*)a,coeff(i)
         if(i.eq.18)xmu(i-17)=coeff(i)
         if((i.gt.25).and.(i.lt.36))xmu(i-24)=coeff(i)
 30   continue
      close(3)
c
c     ************************************************************
c     background counts/sec needed to be subtracted out from the
c     3 channels
c     ************************************************************
c
      open(4,file='background.dat',status='old')
      read(4,*)bkg1,bkg2,bkg3
      close(4)
c
c     ************************************************************
c     The coefficients of the carbon degradation curve. The curve
c     is carbon deposit = acoeff + bcoeff*exp(-day/ccoeff) in (A).
c     ************************************************************
c
      open(10,file='carbon.dat',status='old')
      read(10,*)acoeff,bcoeff,ccoeff
      close(10)
c
c     ************************************************************
c     The beta's, beta1(i) and beta3(i), for channels 1 and 3 are 
c     calculated as a first step to calculate the flux. 
c     ************************************************************
c
      sum1=0.0
      sum3=0.0
      do 13 i=1,11
         if(i.ne.2)sum1=sum1+sol(i)
         sum3=sum3+sol(i)
 13   continue
c
      do 23 i=1,11
         beta1(i)=0.0
         beta3(i)=0.0
         if(i.ne.2)beta1(i)=sol(i)/sum1
         beta3(i)=sol(i)/sum3
 23   continue
c
c     ************************************************************
c     calculate the sums 260 to 340 and 100 to 500
c     ************************************************************
c
      s27=0.0
      do 22 i=27,34
         s27=s27+flux(i)
 22   continue
      s11=0.0
      do 24 i=11,itot
         s11=s11+flux(i)
 24   continue
c
c     ************************************************************
c     read the first channel and second channel count rates
c     ************************************************************
c
      write(*,*)'filename='
      read(*,'(a)')filename
      ii=1
      open(14,file=filename,status='old')
 55   read(14,888,end=901)iyr,iday(ii),jsec(ii),ch1a(ii),sd1(ii),
     +    ch2(ii),sd2(ii),ch1b(ii),sd3(ii),x1(ii),y1(ii),z1(ii),
     +    r1(ii),r(ii),fold,gold
c
      time(ii)=iday(ii)-1+jsec(ii)/86400.0
c
      if(iyr.eq.1996)totday=0
      if(iyr.eq.1997)totday=366
      if(iyr.eq.1998)totday=366+365
      if(iyr.eq.1999)totday=366+2*365
      if(iyr.eq.2000)totday=366+3*365
      if(iyr.eq.2001)totday=366+3*365+366
      if(iyr.eq.2002)totday=366+3*365+366+365
      if(iyr.eq.2003)totday=366+3*365+366+2*365
      if(iyr.eq.2004)totday=366+3*365+366+3*365
      if(iyr.eq.2005)totday=366+3*365+366+3*365+366
      if(iyr.eq.2006)totday=366+3*365+366+3*365+366+365
      if(iyr.eq.2007)totday=366+3*365+366+3*365+366+2*365
      if(iyr.eq.2008)totday=366+3*365+366+3*365+366+3*365
      if(iyr.eq.2009)totday=366+3*365+366+3*365+366+3*365+366
      if(iyr.eq.2010)totday=366+3*365+366+3*365+366+3*365+366+365
c
      day=time(ii)+16.0+totday
      arg=-day/ccoeff
      j=ii
c     carbon deposit is calculated here.
      carb(j)=acoeff+bcoeff*exp(arg)
      ii=ii+1
      goto 55
 901  continue
      close(14)
      isave=ii-1
 888  format(1x,i4,1x,i3,1x,i5,1x,f8.3,1x,e10.3,1x,f10.3,
     +1x,e10.3,1x,f8.3,1x,e10.3,4(1x,e13.5),1x,f8.6,2(e13.5))
c
c     ***********************************************************
c     The output file is named as semflux.dat
c     ***********************************************************
c
      open(50,file='semflux.dat',status='new')
c
      do 100 i=1,isave
c
         if(ch1a(i).eq.0.0)then
            fl1(i)=0.0
            fl2=0.0
            goto 797
         end if
c
c     ************************************************************
c     calculate the normalising factor and the averaged first 
c     order 26-34 num solar flux.
c     ************************************************************
c
         sum1=0.0
         sum3=0.0
         do 110 lk=1,11
            arg1=xmu(lk)*carb(i)
            arg3=xmu(lk)*carb(i)
            trans1=exp(arg1)
            trans3=exp(arg3)
            sum1=sum1+beta1(lk)*q1(lk)*trans1
            sum3=sum3+beta3(lk)*q3(lk)*trans3
 110     continue
         fl1old(i)=((ch1a(i)-bkg1)*r(i)*r(i)/sum1)/area
         fl3old(i)=((ch1b(i)-bkg3)*r(i)*r(i)/sum3)/area
         fl1new(i)=0.8818843565*fl1old(i)
         fl3new(i)=0.814055698*fl3old(i)
         fl1(i)=(0.891*(ch1a(i)-bkg1)*r(i)*r(i)/sum1
     +          +0.814*(ch1b(i)-bkg3)*r(i)*r(i)/sum3)/2.0/area
         fac=fl1(i)/s27
c
         if(eff(25).eq.0)then
            fl2=0.0
            goto 797
         end if
c
c
c     *************************************************************
c     normalize the solar flux longwards of 100 A
c     *************************************************************
c
         do 95 ik=11,itot
            fluxn(ik)=flux(ik)*fac
 95      continue
c
c     *************************************************************
c     calculate the transmission factor
c     *************************************************************
c
         do 90 ik=1,itot
            arg=coeff(ik)*carb(i)
            trans(ik)=exp(arg)
 90      continue
c
c     *************************************************************
c     calculate alpha which measures the deviation of the 0.1 to 10
c     nm solar flux from the composite flux spectrum.
c     *************************************************************
c
         sum11=0.0
         do 85 ik=11,itot
            sum11=sum11+fluxn(ik)*eff(ik)*trans(ik)
 85      continue
         sum0=0.0
         do 80 ik=1,10
            sum0=sum0+flux(ik)*eff(ik)*trans(ik)
 80      continue
         alpha=(((ch2(i)-bkg2)*r(i)*r(i)/area) - sum11)/sum0
c
c     *************************************************************
c     calculate the channel 2 xk coefficient
c     *************************************************************
c
         tot=0.0
         do 70 ik=1,itot
            if(ik.gt.10)fl=fluxn(ik)
            if(ik.lt.11)fl=alpha*flux(ik)
            tot=tot+fl
 70      continue
         sumk=0.0
         do 75 ik=1,itot
            if(ik.gt.10)fl=fluxn(ik)
            if(ik.lt.11)fl=alpha*flux(ik)
            sumk=sumk+fl*eff(ik)*trans(ik)/tot
 75      continue
            xk=1/(area*sumk)
c
c     *************************************************************
c     channel 2 flux, fl2, is calculated
c     *************************************************************
c
         fl2=xk*(ch2(i)-bkg2)*r(i)*r(i)
c
         fl3=fl2 - s11*fac
         fl4=fl2-fl3
         if(fl1(i).lt.0)goto 799
         if(fl3.lt.0)goto 799
 797     write(50,800)iyr,iday(i),jsec(i),ch1a(i),sd1(i),
     +    ch2(i),sd2(i),ch1b(i),sd3(i),x1(i),y1(i),z1(i),
     +    r1(i),r(i),fl1(i),fl2

 799     crap=0.0
 100  continue
c      
 800  format(1x,i4,1x,i3,1x,i5,1x,f8.3,1x,e10.3,1x,f10.3,
     +1x,e10.3,1x,f8.3,1x,e10.3,4(1x,e13.5),1x,f8.6,e13.5,e13.5)
       stop
       end
