;.r Venus_combine_data.pro
;
;
;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
;
; **************************** SECTION 1 - GAUSSIAN FUNCITONS ***********************************************
;
; Gaussian functions defined and called on in the code for fitting. 
; Two functions are calculated, one for a blended Gaussian
; consisting of two Gaussians, and another of a single Gaussian 
;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


;------------------------------------------------------------------------------------------
; function gfunct
; define Gaussian function for deblend case.
; This will be called later on in the code and the variable F will be
; returned.
; A[] array is defined near the bottom of the code in:
;    parinfo[*].value = [2.31e-13, 5577.01, 0.079, 2.11e-12, 5577.232];    **** UPDATE VALUES *****
;
; A[0] = Amplitude of Venusian green line
; A[1] = central wavelength of of Venusian green line
; A[2] = sigma (spread in the data, defined by the dispersion of the instrument, same for both lines)
; A[3] = Amplitude of the Terrestrial green line
; A[4] = central wavelength of the Terrestrial green line
; 
;------------------------------------------------------------------------------------------


    function gfunct,lambda, A, F, pder ;define Gaussian function for deblend case

    F=A[0]*exp(-(lambda-A[1])^2.0/(2.0*(A[2]^2.0)))+A[3]*exp(-(lambda-A[4])^2.0/(2.0*(A[2]^2.0)))

    IF N_PARAMS() GE 4 THEN $
        pder = [[exp(-(lambda-A[1])^2/(2*(A[2]^2)))], [(lambda-A[1])*A[0]/A[2]^2*exp(-(lambda-A[1])^2/(2*(A[2]^2)))], $
                [(lambda-A[1])*A[0]/A[2]^3*exp(-(lambda-A[1])^2/(2*(A[2]^2)))]]

    return, F
END


;------------------------------------------------------------------------------------------
; function gfunct2
; define Gaussian function for individual line
;------------------------------------------------------------------------------------------

    function gfunct2,lambda, A, F, pder ;define Gaussian function for individual line

    F=A[0]*exp(-(lambda-A[1])^2.0/(2.0*(A[2]^2.0)))

    IF N_PARAMS() GE 4 THEN $
        pder = [[exp(-(lambda-A[1])^2/(2*(A[2]^2)))], [(lambda-A[1])*A[0]/A[2]^2*exp(-(lambda-A[1])^2/(2*(A[2]^2)))],$
                [(lambda-A[1])*A[0]/A[2]^3*exp(-(lambda-A[1])^2/(2*(A[2]^2)))]]
    return, F
    END

;------------------------------------------------------------------------------------------
; function OIflux
; A[0]=CO2 flux contribution, A[1]=co2red, A[2]=h2ored, A[3]=co2green, A[4]=h2ogreen
;------------------------------------------------------------------------------------------

    function OIflux, lineratio, A, F; A[0]=CO2 flux contribution, A[1]=co2red, A[2]=h2ored, A[3]=co2green, A[4]=h2ogreen
    B = A[1] / (A[2]+A[1])
    C = (lineratio*A[2] - A[4]) / (A[3] - lineratio*A[1])
    F = A[0]/(B*C)
    return, F 
    END


;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
;
; ****************************** SECTION 2 - READING IN DATA  ***********************************************

;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

READCOL, 'combine_images_inputs.dat', index, position, month, day, year, airmass, exposure, relative_vel, FORMAT='A,A,A,A,I,F,F,F'



;-------------------
;Equator exposures
;-------------------
READCOL,'output_0012_flux.dat', lambda1_flux, flux1_flux, flux1_flux_err, FORMAT='F,F,F'
READCOL,'output_0012_rayleigh.dat', lambda1_rayleigh, flux1_rayleigh, flux1_rayleigh_err, FORMAT='F,F,F'
READCOL,'output_0012_normalized.dat', lambda1_norm, flux1_norm, flux1_norm_err, FORMAT='F,F,F'

READCOL,'output_0013_flux.dat', lambda2_flux, flux2_flux, flux2_flux_err, FORMAT='F,F,F'
READCOL,'output_0013_rayleigh.dat', lambda2_rayleigh, flux2_rayleigh, flux2_rayleigh_err, FORMAT='F,F,F'
READCOL,'output_0013_normalized.dat', lambda2_norm, flux2_norm, flux2_norm_err, FORMAT='F,F,F'



;-------------------
;South exposures
;-------------------
READCOL,'output_0006_flux.dat', lambda3_flux, flux3_flux, flux3_flux_err, FORMAT='F,F<F'
READCOL,'output_0006_rayleigh.dat', lambda3_rayleigh, flux3_rayleigh, flux3_rayleigh_err, FORMAT='F,F,F'
READCOL,'output_0006_normalized.dat', lambda3_norm, flux3_norm, flux3_norm_err, FORMAT='F,F,F'

READCOL,'output_0011_flux.dat', lambda4_flux, flux4_flux, flux4_flux_err, FORMAT='F,F,F'
READCOL,'output_0011_rayleigh.dat', lambda4_rayleigh, flux4_rayleigh, flux4_rayleigh_err, FORMAT='F,F,F'
READCOL,'output_0011_normalized.dat', lambda4_norm, flux4_norm, flux4_norm_err, FORMAT='F,F,F'



lambda = lambda1_flux




;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
;
; ****************************** SECTION 3 - Combine Data  **************************************************

;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


;-----------------------------------------------------------------------------------------------------------
; Combine a straight average of the data, no weighting, as they should
; all have the same exposure time.
;-----------------------------------------------------------------------------------------------------------


equator_flux =     (flux1_flux + flux2_flux)/2.0
equator_rayleigh = (flux1_rayleigh + flux2_rayleigh)/2.0
;equator_norm =     (flux1_norm + flux2_norm)/2.0
equator_norm =     (flux1_norm/MAX(flux1_norm) + flux2_norm/MAX(flux2_norm))/2.0


south_flux =     (flux3_flux + flux4_flux)/2.0
south_rayleigh = (flux3_rayleigh + flux4_rayleigh)/2.0
;south_norm =     (flux3_norm + flux4_norm)/2.0
south_norm =     (flux2_norm/MAX(flux3_norm) + flux4_norm/MAX(flux4_norm))/2.0


;-----------------------------------------------------------------------------------------------------------
; Calculate error for combined spectra, using the propagation:
; error = sqrt[ (flux1_err)^2 + (flux2_err)^2 + ...(fluxn_err)^2 ]
;;-----------------------------------------------------------------------------------------------------------

equator_flux_err     = sqrt((flux1_flux_err)^2.0     + (flux2_flux_err)^2.0)
equator_rayleigh_err = sqrt((flux1_rayleigh_err)^2.0 + (flux2_rayleigh_err)^2.0)
equator_norm_err     = sqrt((flux1_norm_err)^2.0     + (flux2_norm_err)^2.0)

south_flux_err     = sqrt((flux1_flux_err)^2.0     + (flux2_flux_err)^2.0)
south_rayleigh_err = sqrt((flux1_rayleigh_err)^2.0 + (flux2_rayleigh_err)^2.0)
south_norm_err     = sqrt((flux1_norm_err)^2.0     + (flux2_norm_err)^2.0)



;-----------------------------------------------------------------------------------------------------------
; Calculate the Doppler shift and set a range for the Venusian green line
; as well an inital guess for the position
; Calculate max and min of whole data set
;-----------------------------------------------------------------------------------------------------------
 
  doppler_shift = relative_vel[0] * 5577.3 / 3e5


   Venus_lambda_guess = doppler_shift + 5577.3   ; Set Veusian greenline wavelength guess
   Venus_lambda_upper = 5577.3 + doppler_shift + 0.08  ; Upper wavelength bracket around emission
   Venus_lambda_lower = 5577.3 + doppler_shift - 0.08  ; Lower wavelegnth braket around emission


;-----------------------------------------------------------------------------------------------------------
; Calculate max and min of whole data set
;-----------------------------------------------------------------------------------------------------------
   equator_flux_max         = MAX(equator_flux, index)
   equator_flux_max_err     = equator_flux_err(index)

   equator_rayleigh_max     = MAX(equator_rayleigh, index)
   equator_rayleigh_max_err = equator_rayleigh_err(index)

   equator_norm_max         = MAX(equator_norm, index)
   equator_norm_max_err     = equator_norm_err(index)


   equator_flux_min         = MIN(equator_flux, index)
   equator_flux_min_err     = equator_flux_err(index)

   equator_rayleigh_min     = MIN(equator_rayleigh, index)
   equator_rayleigh_min_err = equator_rayleigh_err(index)

   equator_norm_min         = MIN(equator_norm, index)
   equator_norm_min_err     = equator_norm_err(index)



   south_flux_max         = MAX(south_flux, index)
   south_flux_max_err     = south_flux_err(index)

   south_rayleigh_max     = MAX(south_rayleigh, index)
   south_rayleigh_max_err = south_rayleigh_err(index)

   south_norm_max         = MAX(south_norm, index)
   south_norm_max_err     = south_norm_err(index)


   south_flux_min         = MIN(south_flux, index)
   south_flux_min_err     = south_flux_err(index)

   south_rayleigh_min     = MIN(south_rayleigh, index)
   south_rayleigh_min_err = south_rayleigh_err(index)

   south_norm_min         = MIN(south_norm, index)
   south_norm_min_err     = south_norm_err(index)



;-----------------------------------------------------------------------------------------------------------
; Find the max value for the terrestrial and the venusian peaks in
; flux, rayleigh, and normalized units and their wavelength
;-----------------------------------------------------------------------------------------------------------

   linerange    = WHERE((lambda GE 5577.23)AND(lambda LE 5577.4))   ;Specify wavelength range we will analyze

   equator_terr_flux_max     = MAX(equator_flux(linerange),index)
   equator_terr_flux_max_err =     equator_flux_err(linerange(index))
   equator_terr_lambda_max   =     lambda1_flux(linerange(index))

   equator_terr_rayleigh_max = MAX(equator_rayleigh(linerange),index)
   equator_terr_rayleigh_err =     equator_rayleigh_err(linerange(index))

   equator_terr_norm_max     = MAX(equator_norm(linerange),index)
   equator_terr_norm_max_err =     equator_norm_err(linerange(index))


   south_terr_flux_max     = MAX(south_flux(linerange),index)
   south_terr_flux_max_err =     south_flux_err(linerange(index))
   south_terr_lambda_max   =     lambda1_flux(linerange(index))

   south_terr_rayleigh_max = MAX(south_rayleigh(linerange),index)
   south_terr_rayleigh_err =     south_rayleigh_err(linerange(index))

   south_terr_norm_max     = MAX(south_norm(linerange),index)
   south_terr_norm_max_err =     south_norm_err(linerange(index))



   linerange    = WHERE((lambda GE Venus_lambda_lower) AND (lambda LE Venus_lambda_upper))   ;Specify wavelength range we will analyze

   equator_venus_flux_max     = MAX(equator_flux(linerange),index)
   equator_venus_flux_max_err =     equator_flux_err(linerange(index))
   equator_venus_lambda_max   =     lambda1_flux(linerange(index))

   equator_venus_rayleigh_max = MAX(equator_rayleigh(linerange),index)
   equator_venus_rayleigh_err =     equator_rayleigh_err(linerange(index))

   equator_venus_norm_max     = MAX(equator_norm(linerange),index)
   equator_venus_norm_max_err =     equator_norm_err(linerange(index))


   south_venus_flux_max     = MAX(south_flux(linerange),index)
   south_venus_flux_max_err =     south_flux_err(linerange(index))
   south_venus_lambda_max   =     lambda1_flux(linerange(index))

   south_venus_rayleigh_max = MAX(south_rayleigh(linerange),index)
   south_venus_rayleigh_err =     south_rayleigh_err(linerange(index))

   south_venus_norm_max     = MAX(south_norm(linerange),index)
   south_venus_norm_max_err =     south_norm_err(linerange(index))


;-----------------------------------------------------------------------------------------------------------
; Set a wavelength for all the points to be plotte again
; Calculate the average airmass for combined equator and south limb exposure.
;-----------------------------------------------------------------------------------------------------------

lambda =  lambda1_flux

airmass_equator = (airmass[0]+airmass[1])/2.0 ;average airmass of combined exposures
airmass_south   = (airmass[2]+airmass[3])/2.0 ;will need for accounting for terrestrial Rayleigh scaled to 1 airmass





;------------------------------------------------------------------------------------------------------------------------------------
;------------------------------------------------------------------------------------------------------------------------------------
;------------------------------------------------------------------------------------------------------------------------------------


   ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   ;
   ; ****************************************** FLUX UNITS EQUATOR *****************************************
   ;
   ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   PRINT,''
   PRINT, ''

;   tmpread=''
;   READ,'Enter 1 for Flux units, 2 for Rayleighs, or 3 for relative scale: ',tmpread


;   If(tmpread EQ '1')Then Begin

        PRINT,''
        PRINT, ''
        print, 'The following calculations will be done in cgs flux units for the equator'
        PRINT,''
        PRINT, ''

;        Hartley_Oct1     = [equator_flux]

       ;------------------------------------------------------------------------------------------
       ; define parameters to be fit THIS WAS CHANGED A LOT!!!!
       ; in the last parinfo, set the guessed flux, Venusian greenline
       ; position, dispersion of the instrument at that wavelength (in Angstroms), guessed
       ; telluric flux, position, and dispersion at that wavelength
       ;
       ; In the updated version of this code that Kenza altered, parinfo was
       ; changed from fixed values to "limits" and an initial guess was
       ; made.  This means that for each parameter (there are 5, the flux of
       ; Venusian, wavelength of Venusian, dispersion, flux of telluric, and
       ; wavelength of telluric), a rough guess is made on each, just to
       ; start somewhere, and then the fit is run across a range which is
       ; defined in you parinfo[#].limits.
       ;
       ; For example, the range of the first parameter (peak flux of Venusian
       ; line) is set to 2.13e-13, and the range is set to 1e-20 through
       ; 1e-10.  The numbers need to be defined as doubles (IDL defaults to
       ; float) or else not enough values past the decimal will be considered.

       ;    parinfo = replicate({value:0.D, fixed:0, limited:[1,1], limits:[0.,1e6]},5)
       ;    parinfo = replicate({value:0.D, fixed:1, limited:[0,0], limits:[0.0,0.0]},6)
       ;    parinfo[0].fixed = 0  Venusian guessed flux
       ;    parinfo[1].fixed = 0  Venusian line position
       ;    parinfo[3].fixed = 0  Dispersion at that order (0.0709 for 47th  order of echelle)
       ;    parinfo[3].fixed = 0  Terrestrial wavelength
       ;    parinfo[4].fixed = 0  Terrestrial flux
       ;    parinfo[*].value  = [2.31e-13, 5577.012, 0.079, 2.11e-12, 5577.232];          

       ;    parinfo[*].value = [4.5e-13, 5577.01, 0.0709, 2.3e-12, 5577.23, 0.0709];   
       ;    parinfo[*].value = [2e-14, 5577.13, 0.0709, 9e-14, 5577.31, 0.0709];       

       ;------------------------------------------------------------------------------------------

       ;Set limits for your parameters and give an initial guess

         parinfo = replicate({fixed:0, limited:[1,1], limits:[1e-20,1e6]},5)
         parinfo[0].limits = [1e-20,1e-08]     ; Venusian Flux
;         parinfo[1].limits = [5577.4, 5577.58]
         parinfo[1].limits = [DOUBLE(Venus_lambda_lower[0]), DOUBLE(Venus_lambda_upper[0])]
         parinfo[2].limits = [0.078,0.08]
         parinfo[3].limits = [1e-20,1e-8]     ; Terrestrial Flux
         parinfo[4].limits = [5577.23, 5577.4]
     
         start             = [DOUBLE(equator_venus_flux_max), DOUBLE(Venus_lambda_guess[0]), 0.079D, DOUBLE(equator_terr_flux_max), 5577.3D] ;INITIAL GUESSES


         P=fltarr(5,1);5 is number of parameters, 1 is number of spectra to be analyzed
         perror=fltarr(5,1)
    
      ;-------------------------------------------------------------------------------------------
      ; this used to have 2048 instead of N_ELEMENTS but since the array was
      ; changed to just be the wavelength range of the Doppler shifted
      ; Venusian and telluric green line, and I may slightly alter that range
      ; values for other data sets, the number of elements was just written
      ; as the same number of Venus_solsub (remember  that the number of
      ; element for Venus_solsub was changed right before this section
      ;-------------------------------------------------------------------------------------------

        fluxfit     = fltarr(N_ELEMENTS(equator_flux),1)  ;this used to have 2048 instead of N_ELEMENTS
        fluxfittell = fltarr(N_ELEMENTS(equator_flux),1)
        residual    = fltarr(N_ELEMENTS(equator_flux),1)
        totresidual = fltarr(1)
        chisquared  = fltarr(1)
        redchi      = fltarr(1)
        averesid    = fltarr(1)
        ew5577      = fltarr(1)
        ewerr5577   = fltarr(1)


      ;------------------------------------------------------------------------------------------------------------
      ; ------------------------------------ PLOTTING IN FLUX UNITS EQUATOR----------------------------------------
      ;------------------------------------------------------------------------------------------------------------


        angstrom = '!6!sA!r!u!9 %!6!n' ;define the angstrom character

        print, 'begin loop'

        LOADCT,39,/SILENT
        !P.CHARTHICK=2.5
        !P.POSITION=[.17,.12,.87,.95] ;set lower left and upper right corners of plotting box
        !X.CHARSIZE = 1.3 
        !Y.CHARSIZE = 1.3
        !X.THICK=4
        !Y.THICK=4

        SET_PLOT,'PS'
  ;      DEVICE,FILENAME='16July2012_greenline_equator_combined_flux_FINAL.ps',/COLOR,/LANDSCAPE,YSIZE=7,XSIZE=9.5,/COURIER,/INCHES
        DEVICE,FILENAME=strcompress(string(day[0])+month[0]+'2012_greenline_'+position[0]+'_combined_flux_FINAL.ps',/remove_all),$
          /COLOR,/LANDSCAPE,YSIZE=7,XSIZE=9.5,/COURIER,/INCHES


      ;-------------------------------------------------------------------------------------------------------------
      ;--------------------------------------- BEGIN GAUSSIAN FIT FOR PLOT -----------------------------------------
      ;-------------------------------------------------------------------------------------------------------------

;        for i=0, 0 do begin ;0-9 is 10, the number of spectra that are considered. For one spectra, loop goes from 0 to 0

 ;       P[*,i] = MPFITFUN('gfunct', lambda, Hartley_Oct1[*,i], Hartley_Oct1_err[*,i], start, PARINFO=parinfo, $
 ;          NITER=niter, PERROR=perror, STATUS=status, BESTNORM=bestnorm, /NAN) ;, quiet=1)

        P[*] = MPFITFUN('gfunct', lambda, equator_flux[*], equator_flux_err[*], start, PARINFO=parinfo, NITER=niter, PERROR=perror, STATUS=status, BESTNORM=bestnorm, /NAN) ;, quiet=1)

           plot, lambda1_flux, equator_flux[*], xrange=[5576.5, 5578.35], /xstyle, $
                  yrange=[equator_flux_min - equator_flux_min_err*1.2, equator_flux_max + equator_flux_max_err*1.2], $
                  ystyle=1, psym=-4,xtitle='Wavelength ('+angstrom+' )', $ 
                  ytitle='Flux (ergs cm!u-2 !n s!u-1 !n '+angstrom+'!u-1!n)', title=strcompress(month[0]+' '+string(day[0])+', 2012 - '+position[0]+' Limb'), $
                  charsize=1.5, THICK=8
        
    
            ;---------------------------------------------------------------------------------------
            ; Plot error on each point.
            ; If the errors are off the scale, multiple by a small number 
            ; just to compare relative sizes of errors on each point.
            ;--------------------------------------------------------------------------------------
 



             oploterror,lambda1_flux, equator_flux[*], equator_flux_err[*], thick = 4 ;plots the error.    
              wavelength = findgen(3000)/1000 + 5576

 
            ;---------------------------------------------------------------------------------------
            ; Plots Guassian of terrestrail green line:
            ;    P[0,i] * exp(-(wavelength-P[1,i])^2.0 / (2.0 * (P[2,i]^2.0))
            ;
            ; and the Doppler shifted green line:
            ;    P[3,i] * exp(-(wavelength-P[4,i])^2.0 / (2.0 * (P[5,i]^2.0)))
            ; 
            ; added together.  It's drawn as a dashed line.  If it's
            ; not matched in flux to the observed emission, something went 
            ; wrong somewhere.
            ;-----------------------------------------------------------------------------------------

              oplot, wavelength, ((P[0] * exp(-(wavelength-P[1])^2.0 / (2.0 * (P[2]^2.0))) + $
                                   P[3] * exp(-(wavelength-P[4])^2.0 / (2.0 * (P[2]^2.0))))), $
                                   linestyle=2, COLOR=83,THICK=8
  
   
            ;---------------------------------------------------------------------------------------
            ; Guassian plot of just the terrestral green line, so you can check
            ; how it compares to the observed line. If this line fits better
            ; than the blended line, then you don't have green line emission.
            ;---------------------------------------------------------------------------------------
       
              oplot, wavelength, (P[3]*exp(-(wavelength-P[4])^2.0/(2.0*(P[2]^2.0)))), $ 
                                 linestyle=4,COLOR=250,THICK=8 ;Terrestrail Gaussian only

              print, 'parameters:', P
              print, 'parameter errors:', perror


            ;------------------------------------------------------------------------------------------
            ; Integrate over fit and calculate errors
            ;
            ; fluxfit - fits the gaussian of the Venusian line flux.  
            ; This will be zero everywhere except near the emission.  
            ; The max value it should have is that of the guessed emission
            ; strength, which is written in     
            ;    parinfo[*].value = [8e-14, 5577.07, 0.0709, 1e-12, 5577.31, 0.0709]; 
            ;
            ; fluxfittell  - fits Gaussian flux of telluric green line flux
            ; Same hold true here for the zero values everywhere except 
            ; around the line.
            ;
            ; "The residual printed is in units of the 1-sigma error bars, so the residual
            ; should be around 1 or less for most points."
            ;------------------------------------------------------------------------------------------

              print, ""
              print, "start flux fit"
              print, ""


              fluxfit[*]     = P[0] * exp(-(lambda - P[1])^2.0 / (2.0 * (P[2]^2.0)))
              fluxfittell[*] = P[3] * exp(-(lambda - P[4])^2.0 / (2.0 * (P[2]^2.0)))
              residual[*]    = (fluxfit[*] + fluxfittell[*] - (equator_flux[*])) / (equator_flux_err[*]) 

              fiterr=replicate(0.0,N_ELEMENTS(flux1_flux)) ;!!!
              chisquared  = 0.0
              totresidual = 0.0

              for j=0,N_ELEMENTS(equator_flux)-1 do begin ;range of point in lambda to match greenline region
                  totresidual = totresidual + abs(residual[j])
                  chisquared  = chisquared  + residual[j]^2.0
              endfor

              redchi = (chisquared)/4.0
              ew5577 = 2.50663 * P[0] * P[2]

              if sqrt(redchi) gt 1.0 then begin
                  ewerr5577 = 2.50663 * perror * P[2] * sqrt(redchi)
              endif else begin
                  ewerr5577 = 2.50663 * perror * P[2]
              endelse

              averesid=totresidual/N_ELEMENTS(flux1_flux) ;divided by the number of points you're sampling over (N_elements..)

              print, 'average residual:',   averesid
              print, 'reduced chi-squared:',redchi
              print, 'integrated flux:',    ew5577
              print, 'error in flux:',      ewerr5577


;          endfor 

          ;------------------------------------------------------------------------------------------
          ; Print data to file
          ;------------------------------------------------------------------------------------------
  ;          openw, lun, 'output_'+exp_numb[k]+'_flux.dat', /GET_LUN ;**** UPDATE VARIABLE ****
;
;            for i=0, N_ELEMENTS(lambda)-1 do begin
 ;               printf, lun, lambda1_flux, equator_flux
  ;          endfor
   ;         FREE_LUN, lun

         ;------------------------------------------------------------------------------------------

           DEVICE,/CLOSE



;------------------------------------------------------------------------------------------------------------------------------------
;------------------------------------------------------------------------------------------------------------------------------------
;------------------------------------------------------------------------------------------------------------------------------------


   ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   ;
   ; ************************************** RAYLEIGH UNITS EQUATOR****************************************
   ;
   ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   

        PRINT,''
        PRINT, ''
        print, 'The following calculations will be done in Rayleigh units for the equator'
        PRINT,''
        PRINT, ''

        terr_max = -1
        FOR n=16, 19 DO BEGIN
            if equator_rayleigh[n] GT terr_max THEN terr_max = equator_rayleigh[n]
        ENDFOR

        venus_max = -1
        FOR n=19, 23 DO BEGIN
            if equator_rayleigh[n] GT venus_max THEN venus_max = equator_rayleigh[n]
        ENDFOR


       ;------------------------------------------------------------------------------------------
        ; Set limits for your parameters and give an initial guess



         parinfo = replicate({fixed:0, limited:[1,1], limits:[1e-20,1e6]},5)
         parinfo[0].limits = [0, equator_venus_rayleigh_max*2.0]            ; Venusian Relative
         parinfo[1].limits = [DOUBLE(Venus_lambda_lower[0]), DOUBLE(Venus_lambda_upper[0])]
         parinfo[2].limits = [0.078,0.08]
         parinfo[3].limits = [0., equator_terr_rayleigh_max*2.0]            ; Terrestrial Relative
         parinfo[4].limits = [5577.23, 5577.4]
     
         start             = [DOUBLE(equator_venus_rayleigh_max), DOUBLE(Venus_lambda_guess[0]), 0.079D, $
                              DOUBLE(equator_terr_rayleigh_max), 5577.3D]        



         P=fltarr(5,1);5 is number of parameters, 1 is number of spectra to be analyzed
         perror=fltarr(5,1)
    

        fluxfit     = fltarr(N_ELEMENTS(equator_rayleigh),1)  ;this used to have 2048 instead of N_ELEMENTS
        fluxfittell = fltarr(N_ELEMENTS(equator_rayleigh),1)
        residual    = fltarr(N_ELEMENTS(equator_rayleigh),1)
        totresidual = fltarr(1)
        chisquared  = fltarr(1)
        redchi      = fltarr(1)
        averesid    = fltarr(1)
        ew5577      = fltarr(1)
        ewerr5577   = fltarr(1)


      ;------------------------------------------------------------------------------------------------------------
      ; -------------------------------- PLOTTING IN RAYLEIGH UNITS EQUATOR----------------------------------------
      ;------------------------------------------------------------------------------------------------------------


        angstrom = '!6!sA!r!u!9 %!6!n' ;define the angstrom character

        print, 'begin loop'

        LOADCT,39,/SILENT
        !P.CHARTHICK=2.5
        !P.POSITION=[.17,.12,.87,.95] ;set lower left and upper right corners of plotting box
        !X.CHARSIZE = 1.3 
        !Y.CHARSIZE = 1.3
        !X.THICK=4
        !Y.THICK=4

        SET_PLOT,'PS'
       DEVICE,FILENAME=strcompress(string(day[0])+month[0]+'2012_greenline_'+position[0]+'_combined_rayleigh_FINAL.ps',/remove_all),$
          /COLOR,/LANDSCAPE,YSIZE=7,XSIZE=9.5,/COURIER,/INCHES


      ;-------------------------------------------------------------------------------------------------------------
      ;--------------------------------------- BEGIN GAUSSIAN FIT FOR PLOT -----------------------------------------
      ;-------------------------------------------------------------------------------------------------------------

;        for i=0, 0 do begin ;0-9 is 10, the number of spectra that are considered. For one spectra, loop goes from 0 to 0

 ;       P[*,i] = MPFITFUN('gfunct', lambda, Hartley_Oct1[*,i], Hartley_Oct1_err[*,i], start, PARINFO=parinfo, $
 ;          NITER=niter, PERROR=perror, STATUS=status, BESTNORM=bestnorm, /NAN) ;, quiet=1)

        P[*] = MPFITFUN('gfunct', lambda, equator_rayleigh[*], equator_rayleigh_err[*], start, PARINFO=parinfo, $
           NITER=niter, PERROR=perror, STATUS=status, BESTNORM=bestnorm, /NAN) ;, quiet=1)

           plot, lambda, equator_rayleigh[*], xrange=[5576.5, 5578.35], /xstyle, $
                  yrange=[equator_rayleigh_min - equator_rayleigh_min_err*1.2, equator_rayleigh_max + equator_rayleigh_max_err*1.2], $
                  ystyle=9, psym=-4,xtitle='Wavelength ('+angstrom+' )', $ 
                  ytitle='Venusian Intesnsity (R/'+angstrom+')', title=strcompress(month[0]+' '+string(day[0])+', 2012 - '+position[0]+' Limb'), $
                  charsize=1.5, THICK=8
        
            ;-----------------------------------------------------------------------------------------------------
            ; Set up different range on right axis for the Terrestrial intensity. Divide by airmass to account for 
            ; Csc z, as Rayleighs are defined at an airmass of 1, we need to scale to 1.
            ; To get the axis to go down and not up, you need to do a XYOutS command and change the orientation by
            ; 90 degrees the opposite way.
            ;-----------------------------------------------------------------------------------------------------

              axis,yaxis=1,yrange=[(equator_rayleigh_min - equator_rayleigh_min_err*1.2)/airmass_equator, (equator_rayleigh_max + equator_rayleigh_max_err*1.2)/airmass_equator], $
                   ystyle=1, charsize=1.5 
              xloc = !X.Window[1] + 0.09
              yloc = (0.9 - 0.15) / 2.0 + 0.15
              XYOutS, xloc, yloc, 'Terrestrial Intensity (R/' +angstrom+')', Alignment=0.5, Orientation=-90.0, charsize=2.0, /Normal
        


    
            ;---------------------------------------------------------------------------------------
            ; Plot error on each point.
            ; If the errors are off the scale, multiple by a small number 
            ; just to compare relative sizes of errors on each point.
            ;--------------------------------------------------------------------------------------
 
              oploterror,lambda, equator_rayleigh[*], equator_rayleigh_err[*] ;plots the error.    


 
            ;---------------------------------------------------------------------------------------
            ; Plots Guassian of the Doppler shifted Venusian green line:
            ;    P[0,i] * exp(-(wavelength-P[1,i])^2.0 / (2.0 * (P[2,i]^2.0))
            ;
            ; and the terrestrial green line:
            ;    P[3,i] * exp(-(wavelength-P[4,i])^2.0 / (2.0 * (P[5,i]^2.0)))
            ; 
            ; added together.  It's drawn as a dashed line.  If it's
            ; not matched in flux to the observed emission, something went 
            ; wrong somewhere.
            ;-----------------------------------------------------------------------------------------

              wavelength = findgen(3000)/1000 + 5576

              oplot, wavelength, ((P[0] * exp(-(wavelength-P[1])^2.0 / (2.0 * (P[2]^2.0))) + $
                                   P[3] * exp(-(wavelength-P[4])^2.0 / (2.0 * (P[2]^2.0))))), $
                                   linestyle=2, COLOR=83,THICK=8 ; Venusian plus terrestrial
  
   
            ;---------------------------------------------------------------------------------------
            ; Guassian plot of just the terrestral green line, so you can check
            ; how it compares to the observed line. If this line fits better
            ; than the blended line, then you don't have green line emission.
            ;---------------------------------------------------------------------------------------
       
              oplot, wavelength, (P[3]*exp(-(wavelength-P[4])^2.0/(2.0*(P[2]^2.0)))), $ 
                                 linestyle=4,COLOR=250,THICK=8 ;Terrestrail Gaussian only

              
;              oplot, wavelength, (P[0] * exp(-(wavelength-P[1])^2.0 / (2.0 * (P[2]^2.0)))), $
;                                  linestyle=3, COLOR=145,THICK=8 ;Venusian Gaussian only


;             oplot, wavelength, P[3]*exp(-(wavelength-P[4])^2.0/(2.0*(P[2]^2.0))) + P[0] * exp(-(wavelength-P[1])^2.0 / (2.0 * (P[2]^2.0))), linestyle=2, COLOR=83,THICK=8




 

             print, 'parameters:', P
              print, 'parameter errors:', perror


            ;------------------------------------------------------------------------------------------
            ; Integrate over fit and calculate errors
            ;
            ; fluxfit - fits the gaussian of the Venusian line flux.  
            ; This will be zero everywhere except near the emission.  
            ; The max value it should have is that of the guessed emission
            ; strength, which is written in     
            ;    parinfo[*].value = [8e-14, 5577.07, 0.0709, 1e-12, 5577.31, 0.0709]; 
            ;
            ; fluxfittell  - fits Gaussian flux of telluric green line flux
            ; Same hold true here for the zero values everywhere except 
            ; around the line.
            ;
            ; "The residual printed is in units of the 1-sigma error bars, so the residual
            ; should be around 1 or less for most points."
            ;------------------------------------------------------------------------------------------

              print, ""
              print, "start flux fit"
              print, ""


              fluxfit[*]     = P[0] * exp(-(lambda - P[1])^2.0 / (2.0 * (P[2]^2.0)))
              fluxfittell[*] = P[3] * exp(-(lambda - P[4])^2.0 / (2.0 * (P[2]^2.0)))
              residual[*]    = (fluxfit[*] + fluxfittell[*] - (equator_flux[*])) / (equator_flux_err[*]) 

              fiterr=replicate(0.0,N_ELEMENTS(equator_rayleigh)) ;!!!
              chisquared  = 0.0
              totresidual = 0.0

              for j=0,N_ELEMENTS(equator_rayleigh)-1 do begin ;range of point in lambda to match greenline region
                  totresidual = totresidual + abs(residual[j])
                  chisquared  = chisquared  + residual[j]^2.0
              endfor

              redchi = (chisquared)/4.0
              ew5577 = 2.50663 * P[0] * P[2]

              if sqrt(redchi) gt 1.0 then begin
                  ewerr5577 = 2.50663 * perror * P[2] * sqrt(redchi)
              endif else begin
                  ewerr5577 = 2.50663 * perror * P[2]
              endelse

              averesid=totresidual/N_ELEMENTS(equator_rayleigh) ;divided by the number of points you're sampling over (N_elements..)

              print, 'average residual:',   averesid
              print, 'reduced chi-squared:',redchi
              print, 'integrated flux:',    ew5577
              print, 'error in flux:',      ewerr5577


;          endfor 

          ;------------------------------------------------------------------------------------------
          ; Print data to file
          ;------------------------------------------------------------------------------------------
  ;          openw, lun, 'output_'+exp_numb[k]+'_flux.dat', /GET_LUN ;**** UPDATE VARIABLE ****
;
;            for i=0, N_ELEMENTS(lambda)-1 do begin
 ;               printf, lun, lambda1_flux, equator_flux
  ;          endfor
   ;         FREE_LUN, lun

         ;------------------------------------------------------------------------------------------

           DEVICE,/CLOSE








;------------------------------------------------------------------------------------------------------------------------------------
;------------------------------------------------------------------------------------------------------------------------------------
;------------------------------------------------------------------------------------------------------------------------------------

   ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   ;
   ; ************************************** NORMALIZED UNITS EQUATOR****************************************
   ;
   ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   PRINT,''
   PRINT, ''

;   tmpread=''
;   READ,'Enter 1 for Flux units, 2 for Rayleighs, or 3 for relative scale: ',tmpread


;   If(tmpread EQ '1')Then Begin

        PRINT,''
        PRINT, ''
        print, 'The following calculations will be done in normalized units for the equator'
        PRINT,''
        PRINT, ''

        terr_max = -1
        FOR n=10, 15 DO BEGIN
            if equator_norm[n] GT terr_max THEN terr_max = equator_norm[n]
        ENDFOR

        venus_max = -1
        FOR n=15, 20 DO BEGIN
            if equator_norm[n] GT venus_max THEN venus_max = equator_norm[n]
        ENDFOR



       ;------------------------------------------------------------------------------------------
        ; Set limits for your parameters and give an initial guess



         parinfo = replicate({fixed:0, limited:[1,1], limits:[1e-20,1e6]},5)
         parinfo[0].limits = [0, 5]            ; Venusian Relative
         parinfo[1].limits = [DOUBLE(Venus_lambda_lower[0]), DOUBLE(Venus_lambda_upper[0])]
         parinfo[2].limits = [0.078,0.08]
         parinfo[3].limits = [0, 1]            ; Terrestrial Relative
         parinfo[4].limits = [5577.23, 5577.4]
     
         start             = [DOUBLE(equator_venus_norm_max), DOUBLE(Venus_lambda_guess[0]), 0.079D, DOUBLE(equator_terr_norm_max), 5577.29D]        ;INITIAL GUESSES, Relat



         P=fltarr(5,1);5 is number of parameters, 1 is number of spectra to be analyzed
         perror=fltarr(5,1)
    

        fluxfit     = fltarr(N_ELEMENTS(equator_norm),1)  ;this used to have 2048 instead of N_ELEMENTS
        fluxfittell = fltarr(N_ELEMENTS(equator_norm),1)
        residual    = fltarr(N_ELEMENTS(equator_norm),1)
        totresidual = fltarr(1)
        chisquared  = fltarr(1)
        redchi      = fltarr(1)
        averesid    = fltarr(1)
        ew5577      = fltarr(1)
        ewerr5577   = fltarr(1)


      ;------------------------------------------------------------------------------------------------------------
      ; --------------------------------- PLOTTING IN NORMALIZED UNITS EQUATOR-------------------------------------
      ;------------------------------------------------------------------------------------------------------------


        angstrom = '!6!sA!r!u!9 %!6!n' ;define the angstrom character

        print, 'begin loop'

        LOADCT,39,/SILENT
        !P.CHARTHICK=2.5
        !P.POSITION=[.17,.12,.87,.95] ;set lower left and upper right corners of plotting box
        !X.CHARSIZE = 1.3 
        !Y.CHARSIZE = 1.3
        !X.THICK=4
        !Y.THICK=4

        SET_PLOT,'PS'

       DEVICE,FILENAME=strcompress(string(day[0])+month[0]+'2012_greenline_'+position[0]+'_combined_norm_FINAL.ps',/remove_all),$
          /COLOR,/LANDSCAPE,YSIZE=7,XSIZE=9.5,/COURIER,/INCHES


      ;-------------------------------------------------------------------------------------------------------------
      ;--------------------------------------- BEGIN GAUSSIAN FIT FOR PLOT -----------------------------------------
      ;-------------------------------------------------------------------------------------------------------------

;        for i=0, 0 do begin ;0-9 is 10, the number of spectra that are considered. For one spectra, loop goes from 0 to 0

 ;       P[*,i] = MPFITFUN('gfunct', lambda, Hartley_Oct1[*,i], Hartley_Oct1_err[*,i], start, PARINFO=parinfo, $
 ;          NITER=niter, PERROR=perror, STATUS=status, BESTNORM=bestnorm, /NAN) ;, quiet=1)

        P[*] = MPFITFUN('gfunct', lambda, equator_norm[*], equator_norm_err[*], start, PARINFO=parinfo, $
           NITER=niter, PERROR=perror, STATUS=status, BESTNORM=bestnorm, /NAN) ;, quiet=1)

           plot, lambda, equator_norm[*], xrange=[5576.5, 5578.35], /xstyle, $
                  yrange=[equator_norm_min * 1.1, equator_norm_max * 1.1], $
                  ystyle=1, psym=10,xtitle='Wavelength ('+angstrom+' )', $ 
                  ytitle='Normalized Units', title=strcompress(month[0]+' '+string(day[0])+', 2012 - '+position[0]+' Limb'), $
                  charsize=1.5, THICK=8
        
    
            ;---------------------------------------------------------------------------------------
            ; Plot error on each point.
            ; If the errors are off the scale, multiple by a small number 
            ; just to compare relative sizes of errors on each point.
            ;--------------------------------------------------------------------------------------
 
;              oploterror,lambda, equator_norm[*], equator_norm_err[*] ;plots the error.    
              wavelength = findgen(3000)/1000 + 5576

 
            ;---------------------------------------------------------------------------------------
            ; Plots Guassian of terrestrail green line:
            ;    P[0,i] * exp(-(wavelength-P[1,i])^2.0 / (2.0 * (P[2,i]^2.0))
            ;
            ; and the Doppler shifted green line:
            ;    P[3,i] * exp(-(wavelength-P[4,i])^2.0 / (2.0 * (P[5,i]^2.0)))
            ; 
            ; added together.  It's drawn as a dashed line.  If it's
            ; not matched in flux to the observed emission, something went 
            ; wrong somewhere.
            ;-----------------------------------------------------------------------------------------

              oplot, wavelength, ((P[0] * exp(-(wavelength-P[1])^2.0 / (2.0 * (P[2]^2.0))) + $
                                   P[3] * exp(-(wavelength-P[4])^2.0 / (2.0 * (P[2]^2.0))))), $
                                   linestyle=2, COLOR=83,THICK=8
  
   
            ;---------------------------------------------------------------------------------------
            ; Guassian plot of just the terrestral green line, so you can check
            ; how it compares to the observed line. If this line fits better
            ; than the blended line, then you don't have green line emission.
            ;---------------------------------------------------------------------------------------
       
              oplot, wavelength, (P[3]*exp(-(wavelength-P[4])^2.0/(2.0*(P[2]^2.0)))), $ 
                                 linestyle=4,COLOR=250,THICK=8 ;Terrestrail Gaussian only

              print, 'parameters:', P
              print, 'parameter errors:', perror


            ;------------------------------------------------------------------------------------------
            ; Integrate over fit and calculate errors
            ;
            ; fluxfit - fits the gaussian of the Venusian line flux.  
            ; This will be zero everywhere except near the emission.  
            ; The max value it should have is that of the guessed emission
            ; strength, which is written in     
            ;    parinfo[*].value = [8e-14, 5577.07, 0.0709, 1e-12, 5577.31, 0.0709]; 
            ;
            ; fluxfittell  - fits Gaussian flux of telluric green line flux
            ; Same hold true here for the zero values everywhere except 
            ; around the line.
            ;
            ; "The residual printed is in units of the 1-sigma error bars, so the residual
            ; should be around 1 or less for most points."
            ;------------------------------------------------------------------------------------------

              print, ""
              print, "start flux fit"
              print, ""


              fluxfit[*]     = P[0] * exp(-(lambda - P[1])^2.0 / (2.0 * (P[2]^2.0)))
              fluxfittell[*] = P[3] * exp(-(lambda - P[4])^2.0 / (2.0 * (P[2]^2.0)))
              residual[*]    = (fluxfit[*] + fluxfittell[*] - (equator_flux[*])) / (equator_flux_err[*]) 

              fiterr=replicate(0.0,N_ELEMENTS(flux1_norm)) ;!!!
              chisquared  = 0.0
              totresidual = 0.0

              for j=0,N_ELEMENTS(equator_flux)-1 do begin ;range of point in lambda to match greenline region
                  totresidual = totresidual + abs(residual[j])
                  chisquared  = chisquared  + residual[j]^2.0
              endfor

              redchi = (chisquared)/4.0
              ew5577 = 2.50663 * P[0] * P[2]

              if sqrt(redchi) gt 1.0 then begin
                  ewerr5577 = 2.50663 * perror * P[2] * sqrt(redchi)
              endif else begin
                  ewerr5577 = 2.50663 * perror * P[2]
              endelse

              averesid=totresidual/N_ELEMENTS(flux1_norm) ;divided by the number of points you're sampling over (N_elements..)

              print, 'average residual:',   averesid
              print, 'reduced chi-squared:',redchi
              print, 'integrated flux:',    ew5577
              print, 'error in flux:',      ewerr5577


;          endfor 

          ;------------------------------------------------------------------------------------------
          ; Print data to file
          ;------------------------------------------------------------------------------------------
  ;          openw, lun, 'output_'+exp_numb[k]+'_flux.dat', /GET_LUN ;**** UPDATE VARIABLE ****
;
;            for i=0, N_ELEMENTS(lambda)-1 do begin
 ;               printf, lun, lambda1_flux, equator_flux
  ;          endfor
   ;         FREE_LUN, lun

         ;------------------------------------------------------------------------------------------

           DEVICE,/CLOSE



           





 ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 ;
 ;------------------------------------------------------------------------------------------------------------------
 ; ********************************************* SOUTH LIMB PLOTS ***************************************************
 ;------------------------------------------------------------------------------------------------------------------
 ;
 ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



   ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   ;
   ; ****************************************** FLUX UNITS SOUTH*********************************************
   ;
   ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   PRINT,''
   PRINT, ''

;   tmpread=''
;   READ,'Enter 1 for Flux units, 2 for Rayleighs, or 3 for relative scale: ',tmpread


;   If(tmpread EQ '1')Then Begin

        PRINT,''
        PRINT, ''
        print, 'The following calculations will be done in cgs flux units for the south exposures'
        PRINT,''
        PRINT, ''


         parinfo = replicate({fixed:0, limited:[1,1], limits:[1e-20,1e6]},5)
         parinfo[0].limits = [1e-20,1e-10]     ; Venusian Flux
         parinfo[1].limits = [DOUBLE(Venus_lambda_lower[0]), DOUBLE(Venus_lambda_upper[0])]
         parinfo[2].limits = [0.078,0.08]
         parinfo[3].limits = [1e-20,1e-08]     ; Terrestrial Flux
         parinfo[4].limits = [5577.23, 5577.4]


         start             = [DOUBLE(south_venus_flux_max), DOUBLE(Venus_lambda_guess[0]), 0.079D, $
                              DOUBLE(south_terr_flux_max), 5577.3D] ;INITIAL GUESSES
    


         P=fltarr(5,1);5 is number of parameters, 1 is number of spectra to be analyzed
         perror=fltarr(5,1)
    
      ;-------------------------------------------------------------------------------------------
      ; this used to have 2048 instead of N_ELEMENTS but since the array was
      ; changed to just be the wavelength range of the Doppler shifted
      ; Venusian and telluric green line, and I may slightly alter that range
      ; values for other data sets, the number of elements was just written
      ; as the same number of Venus_solsub (remember  that the number of
      ; element for Venus_solsub was changed right before this section
      ;-------------------------------------------------------------------------------------------

        fluxfit     = fltarr(N_ELEMENTS(equator_flux),1)  ;this used to have 2048 instead of N_ELEMENTS
        fluxfittell = fltarr(N_ELEMENTS(equator_flux),1)
        residual    = fltarr(N_ELEMENTS(equator_flux),1)
        totresidual = fltarr(1)
        chisquared  = fltarr(1)
        redchi      = fltarr(1)
        averesid    = fltarr(1)
        ew5577      = fltarr(1)
        ewerr5577   = fltarr(1)


      ;------------------------------------------------------------------------------------------------------------
      ; ------------------------------------------- PLOTTING IN FLUX UNITS ----------------------------------------
      ;------------------------------------------------------------------------------------------------------------


        angstrom = '!6!sA!r!u!9 %!6!n' ;define the angstrom character

        print, 'begin loop'

        LOADCT,39,/SILENT
        !P.CHARTHICK=2.5
        !P.POSITION=[.17,.12,.87,.95] ;set lower left and upper right corners of plotting box
        !X.CHARSIZE = 1.3 
        !Y.CHARSIZE = 1.3   
        !X.THICK=4
        !Y.THICK=4

        SET_PLOT,'PS'

       DEVICE,FILENAME=strcompress(string(day[2])+month[2]+'2012_greenline_'+position[2]+'_combined_flux_FINAL.ps',/remove_all),$
          /COLOR,/LANDSCAPE,YSIZE=7,XSIZE=9.5,/COURIER,/INCHES


      ;-------------------------------------------------------------------------------------------------------------
      ;--------------------------------------- BEGIN GAUSSIAN FIT FOR PLOT -----------------------------------------
      ;-------------------------------------------------------------------------------------------------------------

;        for i=0, 0 do begin ;0-9 is 10, the number of spectra that are considered. For one spectra, loop goes from 0 to 0

 ;       P[*,i] = MPFITFUN('gfunct', lambda, Hartley_Oct1[*,i], Hartley_Oct1_err[*,i], start, PARINFO=parinfo, $
 ;          NITER=niter, PERROR=perror, STATUS=status, BESTNORM=bestnorm, /NAN) ;, quiet=1)

        P[*] = MPFITFUN('gfunct', lambda, south_flux[*], south_flux_err[*], start, PARINFO=parinfo, $
           NITER=niter, PERROR=perror, STATUS=status, BESTNORM=bestnorm, /NAN) ;, quiet=1)

           plot, lambda1_flux, south_flux[*], xrange=[5576.5, 5578.35], /xstyle, $
                 yrange=[south_flux_min - south_flux_min_err*1.2, south_flux_max + south_flux_max_err*1.2], $
                  ystyle=1, psym=-4,xtitle='Wavelength ('+angstrom+' )', $ 
                  ytitle='Flux (ergs cm!u-2 !n s!u-1 !n '+angstrom+'!u-1!n)', title=strcompress(month[0]+' '+string(day[0])+', 2012 - '+position[2]+' Limb'), $
                  charsize=1.5, THICK=8
        
    
            ;---------------------------------------------------------------------------------------
            ; Plot error on each point.
            ; If the errors are off the scale, multiple by a small number 
            ; just to compare relative sizes of errors on each point.
            ;--------------------------------------------------------------------------------------
 
              oploterror,lambda,south_flux[*],south_flux_err[*] ;plots the error.    
              wavelength = findgen(3000)/1000 + 5576

 
            ;---------------------------------------------------------------------------------------
            ; Plots Guassian of terrestrail green line:
            ;    P[0,i] * exp(-(wavelength-P[1,i])^2.0 / (2.0 * (P[2,i]^2.0))
            ;
            ; and the Doppler shifted green line:
            ;    P[3,i] * exp(-(wavelength-P[4,i])^2.0 / (2.0 * (P[5,i]^2.0)))
            ; 
            ; added together.  It's drawn as a dashed line.  If it's
            ; not matched in flux to the observed emission, something went 
            ; wrong somewhere.
            ;-----------------------------------------------------------------------------------------

              oplot, wavelength, ((P[0] * exp(-(wavelength-P[1])^2.0 / (2.0 * (P[2]^2.0))) + $
                                   P[3] * exp(-(wavelength-P[4])^2.0 / (2.0 * (P[2]^2.0))))), $
                                   linestyle=2, COLOR=83,THICK=8
  
   
            ;---------------------------------------------------------------------------------------
            ; Guassian plot of just the terrestral green line, so you can check
            ; how it compares to the observed line. If this line fits better
            ; than the blended line, then you don't have green line emission.
            ;---------------------------------------------------------------------------------------
       
              oplot, wavelength, (P[3]*exp(-(wavelength-P[4])^2.0/(2.0*(P[2]^2.0)))), $ 
                                 linestyle=4,COLOR=250,THICK=8 ;Terrestrail Gaussian only

              print, 'parameters:', P
              print, 'parameter errors:', perror


            ;------------------------------------------------------------------------------------------
            ; Integrate over fit and calculate errors
            ;
            ; fluxfit - fits the gaussian of the Venusian line flux.  
            ; This will be zero everywhere except near the emission.  
            ; The max value it should have is that of the guessed emission
            ; strength, which is written in     
            ;    parinfo[*].value = [8e-14, 5577.07, 0.0709, 1e-12, 5577.31, 0.0709]; 
            ;
            ; fluxfittell  - fits Gaussian flux of telluric green line flux
            ; Same hold true here for the zero values everywhere except 
            ; around the line.
            ;
            ; "The residual printed is in units of the 1-sigma error bars, so the residual
            ; should be around 1 or less for most points."
            ;------------------------------------------------------------------------------------------

              print, ""
              print, "start flux fit"
              print, ""


              fluxfit[*]     = P[0] * exp(-(lambda - P[1])^2.0 / (2.0 * (P[2]^2.0)))
              fluxfittell[*] = P[3] * exp(-(lambda - P[4])^2.0 / (2.0 * (P[2]^2.0)))
              residual[*]    = (fluxfit[*] + fluxfittell[*] - (equator_flux[*])) / (equator_flux_err[*]) 

              fiterr=replicate(0.0,N_ELEMENTS(flux1_flux)) ;!!!
              chisquared  = 0.0
              totresidual = 0.0

              for j=0,N_ELEMENTS(south_flux)-1 do begin ;range of point in lambda to match greenline region
                  totresidual = totresidual + abs(residual[j])
                  chisquared  = chisquared  + residual[j]^2.0
              endfor

              redchi = (chisquared)/4.0
              ew5577 = 2.50663 * P[0] * P[2]

              if sqrt(redchi) gt 1.0 then begin
                  ewerr5577 = 2.50663 * perror * P[2] * sqrt(redchi)
              endif else begin
                  ewerr5577 = 2.50663 * perror * P[2]
              endelse

              averesid=totresidual/N_ELEMENTS(flux1_flux) ;divided by the number of points you're sampling over (N_elements..)

              print, 'average residual:',   averesid
              print, 'reduced chi-squared:',redchi
              print, 'integrated flux:',    ew5577
              print, 'error in flux:',      ewerr5577


;          endfor 

          ;------------------------------------------------------------------------------------------
          ; Print data to file
          ;------------------------------------------------------------------------------------------
  ;          openw, lun, 'output_'+exp_numb[k]+'_flux.dat', /GET_LUN ;**** UPDATE VARIABLE ****
;
;            for i=0, N_ELEMENTS(lambda)-1 do begin
 ;               printf, lun, lambda1_flux, equator_flux
  ;          endfor
   ;         FREE_LUN, lun

         ;------------------------------------------------------------------------------------------

           DEVICE,/CLOSE


;------------------------------------------------------------------------------------------------------------------------------------
;------------------------------------------------------------------------------------------------------------------------------------
;------------------------------------------------------------------------------------------------------------------------------------




  ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   ;
   ; ************************************** RAYLEIGH UNITS SOUTH****************************************
   ;
   ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



        PRINT,''
        PRINT, ''
        print, 'The following calculations will be done in Rayleigh units for the South Limb'
        PRINT,''
        PRINT, ''

        terr_max = -1
        FOR n=16, 19 DO BEGIN
            if equator_rayleigh[n] GT terr_max THEN terr_max = equator_rayleigh[n]
        ENDFOR

        venus_max = -1
        FOR n=19, 23 DO BEGIN
            if equator_rayleigh[n] GT venus_max THEN venus_max = equator_rayleigh[n]
        ENDFOR


       ;------------------------------------------------------------------------------------------
        ; Set limits for your parameters and give an initial guess


        parinfo = replicate({fixed:0, limited:[1,1], limits:[1e-20,1e6]},5)
         parinfo[0].limits = [0, south_venus_rayleigh_max*2.0]            ; Venusian Relative
         parinfo[1].limits = [DOUBLE(Venus_lambda_lower[0]), DOUBLE(Venus_lambda_upper[0])-0.08]
         parinfo[2].limits = [0.078,0.079]
         parinfo[3].limits = [0., south_terr_rayleigh_max*2.0]            ; Terrestrial Relative
         parinfo[4].limits = [5577.23, 5577.4]
     
         start             = [DOUBLE(south_venus_rayleigh_max), DOUBLE(Venus_lambda_guess[0])-0.03, 0.079D, $
                              DOUBLE(south_terr_rayleigh_max), 5577.3D]        



         P=fltarr(5,1);5 is number of parameters, 1 is number of spectra to be analyzed
         perror=fltarr(5,1)
    

        fluxfit     = fltarr(N_ELEMENTS(equator_norm),1)  ;this used to have 2048 instead of N_ELEMENTS
        fluxfittell = fltarr(N_ELEMENTS(equator_norm),1)
        residual    = fltarr(N_ELEMENTS(equator_norm),1)
        totresidual = fltarr(1)
        chisquared  = fltarr(1)
        redchi      = fltarr(1)
        averesid    = fltarr(1)
        ew5577      = fltarr(1)
        ewerr5577   = fltarr(1)


      ;------------------------------------------------------------------------------------------------------------
      ; ------------------------------------------- PLOTTING IN RAYLEIGH UNITS ----------------------------------------
      ;------------------------------------------------------------------------------------------------------------


        angstrom = '!6!sA!r!u!9 %!6!n' ;define the angstrom character

        print, 'begin loop'

        LOADCT,39,/SILENT
        !P.CHARTHICK=2.5
        !P.POSITION=[.17,.12,.87,.95] ;set lower left and upper right corners of plotting box
        !X.CHARSIZE = 1.3 
        !Y.CHARSIZE = 1.3
        !X.THICK=4
        !Y.THICK=4

        SET_PLOT,'PS'

      DEVICE,FILENAME=strcompress(string(day[2])+month[2]+'2012_greenline_'+position[2]+'_combined_rayleigh_FINAL.ps',/remove_all),$
          /COLOR,/LANDSCAPE,YSIZE=7,XSIZE=9.5,/COURIER,/INCHES

      ;-------------------------------------------------------------------------------------------------------------
      ;--------------------------------------- BEGIN GAUSSIAN FIT FOR PLOT -----------------------------------------
      ;-------------------------------------------------------------------------------------------------------------

;        for i=0, 0 do begin ;0-9 is 10, the number of spectra that are considered. For one spectra, loop goes from 0 to 0

 ;       P[*,i] = MPFITFUN('gfunct', lambda, Hartley_Oct1[*,i], Hartley_Oct1_err[*,i], start, PARINFO=parinfo, $
 ;          NITER=niter, PERROR=perror, STATUS=status, BESTNORM=bestnorm, /NAN) ;, quiet=1)

        P[*] = MPFITFUN('gfunct', lambda, south_rayleigh[*], south_rayleigh_err[*], start, PARINFO=parinfo, $
           NITER=niter, PERROR=perror, STATUS=status, BESTNORM=bestnorm, /NAN) ;, quiet=1)

           plot, lambda, south_rayleigh[*], xrange=[5576.5, 5578.35], /xstyle, $
                 yrange=[south_rayleigh_min - south_rayleigh_min_err*1.2, south_rayleigh_max + south_rayleigh_max_err*1.2], $
                  ystyle=9, psym=-4,xtitle='Wavelength ('+angstrom+' )', $ 
                  ytitle='Venusian Intesnsity (R/'+angstrom+')', title=strcompress(month[0]+' '+string(day[0])+', 2012 - '+position[2]+' Limb'), $
                  charsize=1.5, THICK=8
        
            ;-----------------------------------------------------------------------------------------------------
            ; Set up different range on right axis for the Terrestrial intensity. Divide by airmass to account for 
            ; Csc z, as Rayleighs are defined at an airmass of 1, we need to scale to 1.
            ; To get the axis to go down and not up, you need to do a XYOutS command and change the orientation by
            ; 90 degrees the opposite way.
            ;-----------------------------------------------------------------------------------------------------

              axis,yaxis=1,yrange=[(south_rayleigh_min - south_rayleigh_min_err*1.2)/airmass_south, (south_rayleigh_max + south_rayleigh_max_err*1.2)/airmass_south], $
                   ystyle=1, charsize=1.5 
              xloc = !X.Window[1] + 0.09
              yloc = (0.9 - 0.15) / 2.0 + 0.15
              XYOutS, xloc, yloc, 'Terrestrial Intensity (R/' +angstrom+')', Alignment=0.5, Orientation=-90.0, charsize=2.0, /Normal
        


    
            ;---------------------------------------------------------------------------------------
            ; Plot error on each point.
            ; If the errors are off the scale, multiple by a small number 
            ; just to compare relative sizes of errors on each point.
            ;--------------------------------------------------------------------------------------
 
              oploterror,lambda, south_rayleigh[*], south_rayleigh_err[*] ;plots the error.    
              wavelength = findgen(3000)/1000 + 5576

 
            ;---------------------------------------------------------------------------------------
            ; Plots Guassian of the Doppler shifted Venusian green line:
            ;    P[0,i] * exp(-(wavelength-P[1,i])^2.0 / (2.0 * (P[2,i]^2.0))
            ;
            ; and the terrestrial green line:
            ;    P[3,i] * exp(-(wavelength-P[4,i])^2.0 / (2.0 * (P[5,i]^2.0)))
            ; 
            ; added together.  It's drawn as a dashed line.  If it's
            ; not matched in flux to the observed emission, something went 
            ; wrong somewhere.
            ;-----------------------------------------------------------------------------------------

              oplot, wavelength, ((P[0] * exp(-(wavelength-P[1])^2.0 / (2.0 * (P[2]^2.0))) + $
                                   P[3] * exp(-(wavelength-P[4])^2.0 / (2.0 * (P[2]^2.0))))), $
                                   linestyle=2, COLOR=83,THICK=8
  
   
            ;---------------------------------------------------------------------------------------
            ; Guassian plot of just the terrestral green line, so you can check
            ; how it compares to the observed line. If this line fits better
            ; than the blended line, then you don't have green line emission.
            ;---------------------------------------------------------------------------------------
       
              oplot, wavelength, (P[3]*exp(-(wavelength-P[4])^2.0/(2.0*(P[2]^2.0)))), $ 
                                 linestyle=4,COLOR=250,THICK=8 ;Terrestrail Gaussian only

              
;              oplot, wavelength, (P[0] * exp(-(wavelength-P[1])^2.0 / (2.0 * (P[2]^2.0)))), $
;                                  linestyle=3, COLOR=145,THICK=8 ;Venusian Gaussian only


;              oplot, wavelength, P[3]*exp(-(wavelength-P[4])^2.0/(2.0*(P[2]^2.0))) + P[0] * exp(-(wavelength-P[1])^2.0 / (2.0 * (P[2]^2.0))), linestyle=4, COLOR=83,THICK=8




 

             print, 'parameters:', P
              print, 'parameter errors:', perror


            ;------------------------------------------------------------------------------------------
            ; Integrate over fit and calculate errors
            ;
            ; fluxfit - fits the gaussian of the Venusian line flux.  
            ; This will be zero everywhere except near the emission.  
            ; The max value it should have is that of the guessed emission
            ; strength, which is written in     
            ;    parinfo[*].value = [8e-14, 5577.07, 0.0709, 1e-12, 5577.31, 0.0709]; 
            ;
            ; fluxfittell  - fits Gaussian flux of telluric green line flux
            ; Same hold true here for the zero values everywhere except 
            ; around the line.
            ;
            ; "The residual printed is in units of the 1-sigma error bars, so the residual
            ; should be around 1 or less for most points."
            ;------------------------------------------------------------------------------------------

              print, ""
              print, "start flux fit"
              print, ""


              fluxfit[*]     = P[0] * exp(-(lambda - P[1])^2.0 / (2.0 * (P[2]^2.0)))
              fluxfittell[*] = P[3] * exp(-(lambda - P[4])^2.0 / (2.0 * (P[2]^2.0)))
              residual[*]    = (fluxfit[*] + fluxfittell[*] - (equator_flux[*])) / (equator_flux_err[*]) 

              fiterr=replicate(0.0,N_ELEMENTS(south_rayleigh)) ;!!!
              chisquared  = 0.0
              totresidual = 0.0

              for j=0,N_ELEMENTS(south_rayleigh)-1 do begin ;range of point in lambda to match greenline region
                  totresidual = totresidual + abs(residual[j])
                  chisquared  = chisquared  + residual[j]^2.0
              endfor

              redchi = (chisquared)/4.0
              ew5577 = 2.50663 * P[0] * P[2]

              if sqrt(redchi) gt 1.0 then begin
                  ewerr5577 = 2.50663 * perror * P[2] * sqrt(redchi)
              endif else begin
                  ewerr5577 = 2.50663 * perror * P[2]
              endelse

              averesid=totresidual/N_ELEMENTS(south_rayleigh) ;divided by the number of points you're sampling over (N_elements..)

              print, 'average residual:',   averesid
              print, 'reduced chi-squared:',redchi
              print, 'integrated flux:',    ew5577
              print, 'error in flux:',      ewerr5577


;          endfor 

          ;------------------------------------------------------------------------------------------
          ; Print data to file
          ;------------------------------------------------------------------------------------------
  ;          openw, lun, 'output_'+exp_numb[k]+'_flux.dat', /GET_LUN ;**** UPDATE VARIABLE ****
;
;            for i=0, N_ELEMENTS(lambda)-1 do begin
 ;               printf, lun, lambda1_flux, equator_flux
  ;          endfor
   ;         FREE_LUN, lun

         ;------------------------------------------------------------------------------------------

           DEVICE,/CLOSE
















;------------------------------------------------------------------------------------------------------------------------------------
;------------------------------------------------------------------------------------------------------------------------------------
;------------------------------------------------------------------------------------------------------------------------------------




  ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   ;
   ; ************************************** NORMALIZED UNITS SOUTH****************************************
   ;
   ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   PRINT,''
   PRINT, ''

;   tmpread=''
;   READ,'Enter 1 for Flux units, 2 for Rayleighs, or 3 for relative scale: ',tmpread


;   If(tmpread EQ '1')Then Begin

        PRINT,''
        PRINT, ''
        print, 'The following calculations will be done in normalized units for the south limb'
        PRINT,''
        PRINT, ''

;        terr_max = -1
;        FOR n=16, 19 DO BEGIN
;            if south_norm[n] GT terr_max THEN terr_max = south_norm[n]
;        ENDFOR

;        venus_max = -1
;        FOR n=19, 23 DO BEGIN
;            if south_norm[n] GT venus_max THEN venus_max = south_norm[n]
;        ENDFOR



       ;------------------------------------------------------------------------------------------
        ; Set limits for your parameters and give an initial guess


         parinfo = replicate({fixed:0, limited:[1,1], limits:[1e-20,1e6]},5)
         parinfo[0].limits = [0, 5]            ; Venusian Relative
         parinfo[1].limits = [DOUBLE(Venus_lambda_lower[0]), DOUBLE(Venus_lambda_upper[0])]
         parinfo[2].limits = [0.078,0.08]
         parinfo[3].limits = [0, 1.2]            ; Terrestrial Relative
         parinfo[4].limits = [5577.23, 5577.4]
     
         start             = [DOUBLE(south_venus_norm_max),  DOUBLE(Venus_lambda_guess[0]), 0.079D, DOUBLE(south_terr_norm_max), 5577.29D]        ;INITIAL GUESSES, Relat




         P=fltarr(5,1);5 is number of parameters, 1 is number of spectra to be analyzed
         perror=fltarr(5,1)
    

        fluxfit     = fltarr(N_ELEMENTS(south_norm),1)  ;this used to have 2048 instead of N_ELEMENTS
        fluxfittell = fltarr(N_ELEMENTS(south_norm),1)
        residual    = fltarr(N_ELEMENTS(south_norm),1)
        totresidual = fltarr(1)
        chisquared  = fltarr(1)
        redchi      = fltarr(1)
        averesid    = fltarr(1)
        ew5577      = fltarr(1)
        ewerr5577   = fltarr(1)


      ;------------------------------------------------------------------------------------------------------------
      ; ---------------------------------- PLOTTING IN NORMALIZED UNITS SOUTH--------------------------------------
      ;------------------------------------------------------------------------------------------------------------


        angstrom = '!6!sA!r!u!9 %!6!n' ;define the angstrom character

        print, 'begin loop'

        LOADCT,39,/SILENT
        !P.CHARTHICK=2.5
        !P.POSITION=[.17,.12,.87,.95] ;set lower left and upper right corners of plotting box
        !X.CHARSIZE = 1.3 
        !Y.CHARSIZE = 1.3
        !X.THICK=4
        !Y.THICK=4

        SET_PLOT,'PS'
        
      DEVICE,FILENAME=strcompress(string(day[2])+month[2]+'2012_greenline_'+position[2]+'_combined_norm_FINAL.ps',/remove_all),$
          /COLOR,/LANDSCAPE,YSIZE=7,XSIZE=9.5,/COURIER,/INCHES

      ;-------------------------------------------------------------------------------------------------------------
      ;--------------------------------------- BEGIN GAUSSIAN FIT FOR PLOT -----------------------------------------
      ;-------------------------------------------------------------------------------------------------------------

;        for i=0, 0 do begin ;0-9 is 10, the number of spectra that are considered. For one spectra, loop goes from 0 to 0

 ;       P[*,i] = MPFITFUN('gfunct', lambda, Hartley_Oct1[*,i], Hartley_Oct1_err[*,i], start, PARINFO=parinfo, $
 ;          NITER=niter, PERROR=perror, STATUS=status, BESTNORM=bestnorm, /NAN) ;, quiet=1)

        P[*] = MPFITFUN('gfunct', lambda, south_norm[*], south_norm_err[*], start, PARINFO=parinfo, $
           NITER=niter, PERROR=perror, STATUS=status, BESTNORM=bestnorm, /NAN) ;, quiet=1)

           plot, lambda, south_norm[*], xrange=[5576.5, 5578.35], /xstyle, $
                 yrange=[south_norm_min * 1.1, south_norm_max * 1.1], $
                  ystyle=1, psym=10,xtitle='Wavelength ('+angstrom+' )', $ 
                  ytitle='Normalized Units', title=strcompress(month[0]+' '+string(day[0])+', 2012 - '+position[2]+' Limb'), $
                  charsize=1.5, THICK=8
        
    
            ;---------------------------------------------------------------------------------------
            ; Plot error on each point.
            ; If the errors are off the scale, multiple by a small number 
            ; just to compare relative sizes of errors on each point.
            ;--------------------------------------------------------------------------------------
 
;              oploterror,lambda, south_norm[*], south_norm_err[*] ;plots the error.    

              wavelength = findgen(3000)/1000 + 5576

 
            ;---------------------------------------------------------------------------------------
            ; Plots Guassian of terrestrail green line:
            ;    P[0,i] * exp(-(wavelength-P[1,i])^2.0 / (2.0 * (P[2,i]^2.0))
            ;
            ; and the Doppler shifted green line:
            ;    P[3,i] * exp(-(wavelength-P[4,i])^2.0 / (2.0 * (P[5,i]^2.0)))
            ; 
            ; added together.  It's drawn as a dashed line.  If it's
            ; not matched in flux to the observed emission, something went 
            ; wrong somewhere.
            ;-----------------------------------------------------------------------------------------

              oplot, wavelength, ((P[0] * exp(-(wavelength-P[1])^2.0 / (2.0 * (P[2]^2.0))) + $
                                   P[3] * exp(-(wavelength-P[4])^2.0 / (2.0 * (P[2]^2.0))))), $
                                   linestyle=2, COLOR=82,THICK=8
  
   
            ;---------------------------------------------------------------------------------------
            ; Guassian plot of just the terrestral green line, so you can check
            ; how it compares to the observed line. If this line fits better
            ; than the blended line, then you don't have green line emission.
            ;---------------------------------------------------------------------------------------
       
              oplot, wavelength, (P[3]*exp(-(wavelength-P[4])^2.0/(2.0*(P[2]^2.0)))), $ 
                                 linestyle=4,COLOR=250,THICK=8 ;Terrestrail Gaussian only

              print, 'parameters:', P
              print, 'parameter errors:', perror


            ;------------------------------------------------------------------------------------------
            ; Integrate over fit and calculate errors
            ;
            ; fluxfit - fits the gaussian of the Venusian line flux.  
            ; This will be zero everywhere except near the emission.  
            ; The max value it should have is that of the guessed emission
            ; strength, which is written in     
            ;    parinfo[*].value = [8e-14, 5577.07, 0.0709, 1e-12, 5577.31, 0.0709]; 
            ;
            ; fluxfittell  - fits Gaussian flux of telluric green line flux
            ; Same hold true here for the zero values everywhere except 
            ; around the line.
            ;
            ; "The residual printed is in units of the 1-sigma error bars, so the residual
            ; should be around 1 or less for most points."
            ;------------------------------------------------------------------------------------------

              print, ""
              print, "start flux fit"
              print, ""


              fluxfit[*]     = P[0] * exp(-(lambda - P[1])^2.0 / (2.0 * (P[2]^2.0)))
              fluxfittell[*] = P[3] * exp(-(lambda - P[4])^2.0 / (2.0 * (P[2]^2.0)))
              residual[*]    = (fluxfit[*] + fluxfittell[*] - (equator_flux[*])) / (equator_flux_err[*]) 

              fiterr=replicate(0.0,N_ELEMENTS(south_norm)) ;!!!
              chisquared  = 0.0
              totresidual = 0.0

              for j=0,N_ELEMENTS(equator_flux)-1 do begin ;range of point in lambda to match greenline region
                  totresidual = totresidual + abs(residual[j])
                  chisquared  = chisquared  + residual[j]^2.0
              endfor

              redchi = (chisquared)/4.0
              ew5577 = 2.50663 * P[0] * P[2]

              if sqrt(redchi) gt 1.0 then begin
                  ewerr5577 = 2.50663 * perror * P[2] * sqrt(redchi)
              endif else begin
                  ewerr5577 = 2.50663 * perror * P[2]
              endelse

              averesid=totresidual/N_ELEMENTS(Venus_solsub) ;divided by the number of points you're sampling over (N_elements..)

              print, 'average residual:',   averesid
              print, 'reduced chi-squared:',redchi
              print, 'integrated flux:',    ew5577
              print, 'error in flux:',      ewerr5577


;          endfor 

          ;------------------------------------------------------------------------------------------
          ; Print data to file
          ;------------------------------------------------------------------------------------------
  ;          openw, lun, 'output_'+exp_numb[k]+'_flux.dat', /GET_LUN ;**** UPDATE VARIABLE ****
;
;            for i=0, N_ELEMENTS(lambda)-1 do begin
 ;               printf, lun, lambda1_flux, equator_flux
  ;          endfor
   ;         FREE_LUN, lun

         ;------------------------------------------------------------------------------------------

           DEVICE,/CLOSE





STOP
end