;Spectral corrections

FUNCTION get_spec,doy,yr,hours,specheader=specheader
;open a spec file and get the data for specific hours in a day
  IF NOT keyword_set(yr) THEN year=co2_const('cur_year') ELSE year = yr
  year2 = year MOD 100
  yrstr = string(year2,format='(i2.2)')
  doystr = string(doy,format='(i3.3)')
  IF NOT keyword_set(hours) THEN hours = indgen(48) /2.0
  ntimes = n_elements(hours)
  fhead = co2_const('file_head')
  specfile = co2_basedir('spec',yr=year)+fhead+yrstr+'.spec.xdf'
  specs = make_array(5,ntimes)
  specs[0,*] = hours
  specs[1,*] = 1.0
  specs[2,*] = -40.0
  specs[3,*] = 1.0
  specs[4,*] = -40.0
  IF file_test(specfile,/read) THEN BEGIN
    print,'  Reading SPECS from ',specfile
    specfp = assoc_xdf(specfile,header=specheader) ;doy,hhmm,c_cfactor,cflag,q_cfactor,qflag
    specdata = line_xdf(specfp,(doy-1)*48,numlines=48)
    close_xdf,specfp
    goodspec = where((isnotnan(specdata[2,*]) OR isnotnan(specdata[4,*])),ngood,complement=badspec)
    IF ngood GT 0 THEN BEGIN    ;some or all good data - can use
      IF badspec[0] NE -1 THEN BEGIN ;some bad data - median fill
        medc = median(specdata[2,*])
        medq = median(specdata[4,*])
        specdata[2,badspec] = medc
        specdata[4,badspec] = medq
      ENDIF
      timelocs = fix(hours * 2)
      specs[1:4,*] = specdata[2:5,timelocs]
    ENDIF ELSE print,'    No good SPEC data for the day, using 1'
  ENDIF ELSE BEGIN
    print,'  Creating SPEC file ',specfile
    print,'  Using 1 for current specs'
    specarray = make_array(6,17568,/float,value=nan())
    specarray[0,*] = float(fix(findgen(17568) / 48) + 1)
    specarray[1,*] = float(fjday_to_hhmm((findgen(17568) MOD 48) / 48.0))
    specarray[2,*] = 1.0
    specarray[3,*] = -40.0
    specarray[4,*] = 1.0
    specarray[5,*] = -40.0
    specheader = ['doy','hhmm','c_cfactor','cflag','q_cfactor','qflag']
    create_dir,co2_basedir('spec',yr=year)
    write_xdf,specfile,specarray,header=specheader
    specarray = 0b
  ENDELSE
  return,specs
END

FUNCTION calc_spec,dq,dc,dtq,dtc,ndt,dwq,dwc,lagc,noplot=noplot,nospec=nospec,plottitle=plottitle,hz=hz,quiet=quiet
;returns: specs = c_cfactor,cflag,q_cfactor,qflag

;Check to see if arrays are good to use
  IF keyword_set(nospec) THEN BEGIN
    IF NOT keyword_set(quiet) THEN print,'    Skipping Spectral correction (nospec) '
    return,[1,-40,1,-40]
  ENDIF
  lenq=n_elements(dq) & lenc=n_elements(dc) 
  lentq=n_elements(dtq) & lentc=n_elements(dtc)
  lenwq=n_elements(dwq) & lenwc=n_elements(dwc)
  qlen_ok = 0 & clen_ok = 0
  IF (lenq eq lentq) and (lenq eq lenwq) THEN qlen_ok = 1 ELSE BEGIN
    IF NOT keyword_set(quiet) THEN print,'    H2O data length problems'
  ENDELSE
  IF (lenc eq lentc) and (lenc eq lenwc) THEN clen_ok = 1 ELSE BEGIN
    IF NOT keyword_set(quiet) THEN print,'    CO2 data length problems'
  ENDELSE    
  whereq = where(isnan(dq),countq) &  wherec = where(isnan(dc),countc)
  wheretq = where(isnan(dtq),counttq) &  wheretc = where(isnan(dtc),counttc)
  wherewq = where(isnan(dwq),countwq) &  wherewc = where(isnan(dwc),countwc)  
  percentbadq = 100.0 * countq / lenq &  percentbadc = 100.0 * countc / lenc
  percentbadtq = 100.0 * counttq / lentq &  percentbadtc = 100.0 * counttc / lentc
  percentbadwq = 100.0 * countwq / lenwq &  percentbadwc = 100.0 * countwc / lenwc
  q_gooddata = 0 & c_gooddata = 0
  IF (percentbadq LT 20.) AND (percentbadtq LT 20.) AND (percentbadwq LT 20.) THEN q_gooddata=1 ELSE BEGIN
    IF NOT keyword_set(quiet) THEN print,'    More than 20% bad Q vals for fit.'
  ENDELSE
  IF (percentbadc LT 20.) AND (percentbadtc LT 20.) AND (percentbadwc LT 20.) THEN c_gooddata=1 ELSE BEGIN
    IF NOT keyword_set(quiet) THEN print,'    More than 20% bad C vals for fit.'
  ENDELSE
  specs = float([1,((c_gooddata EQ 0)*(-10))+((clen_ok EQ 0)*(-20)),1,((q_gooddata EQ 0)*(-10))+((qlen_ok EQ 0)*(-20))])
  IF n_elements(hz) EQ 0 THEN BEGIN 
    fs = 10.
    hz = 10.
  ENDIF ELSE fs = hz
  badkneeflag=0
  badtspecflag=0

;Finding q_cfactor
  IF (qlen_ok EQ 1) AND (q_gooddata EQ 1) THEN BEGIN
    nq = 2l^long(alog(lenq)/alog(2)) 
    dqc = zapbadval(dq[0:nq-1]) & dtqc = zapbadval(dtq[0:nq-1]) & dwqc = zapbadval(dwq[0:nq-1]) ;fill holes and remove mean
    qmn = mean(dqc,/nan) & tqmn = mean(dtqc,/nan) & wqmn = mean(dwqc,/nan)
    q = dqc - qmn & tq = dtqc - tqmn & wq = dwqc - wqmn
    aq=fft(q,/double)           ;q spectra calcs
    aqhat=double(conj(aq)*aq)
    delfq=(nq-1)/double(fs)
    gqq=2*delfq*aqhat
    atq=fft(tq,/double)         ;tq spectra calcs
    athat=double(conj(atq)*atq)
    gtt=2*delfq*athat
    halfq=nq/2+1                ;extract up to Nyquist freq
    gqq=gqq[0:halfq-2]
    gtt=gtt[0:halfq-2]
    fq=dindgen(halfq-1)*fs/(nq-1)
    gttsm=running_avg(cdr(gtt),11) ;smooth spectra with running average
    gqqsm=running_avg(cdr(gqq),11)
    fsm=cdr(fq)
    gttsmc=gttsm[5:n_elements(gttsm)-6]     ;trim out beginning and end (was halfq-8)
    gqqsmc=gqqsm[5:n_elements(gqqsm)-6]
    fsmc=fsm[5:n_elements(fsm)-6]
    tt=(1/delfq)*total(gttsmc[0:14])  ;normalizing variance (1st 15 pts)
    qq=(1/delfq)*total(gqqsmc[0:14])
    edg=edgedetect(gqqsmc,201)  ;square wavelet to find knee of spectra
    wav=edg[0,*]
    tm=edg[1,*]
    tmind = findi(wav,0.001,~keyword_set(quiet))

    IF isnan(tmind) OR (tmind EQ -1) THEN badkneeflag = 1 ELSE BEGIN 
      IF abs(wav[tmind]-0.001) GT .1 THEN badkneeflag = 1 ELSE IF (tm[tmind]-100) LE 2 THEN badkneeflag = 1
    ENDELSE
    IF badkneeflag EQ 1 THEN BEGIN
      IF NOT keyword_set(quiet) THEN print,'    Hard to find spectra knee. No H2O spectral correction'
      hiind = n_elements(fsmc)  ;for plotting purposes
      hifreq=fsmc[hiind-1] 
   ENDIF ELSE BEGIN
      hiind = tm[tmind] - 100l
      hifreq=fsmc[hiind]
      IF NOT keyword_set(quiet) THEN print,'    Knee found, f = ',hifreq
    ENDELSE
    att=poly_fit(alog10(fsmc),alog10(gttsmc/tt),3) ;fit spectra on whole length
    aqq1flag=0
    aqq=poly_fit(alog10(fsmc[0:hiind-1]),alog10(gqqsmc[0:hiind-1]/qq),2) ;fit upto knee
    IF aqq[2] GT 0.5 THEN BEGIN ;deal with bad fits past knee - convert to first order
      aqq1flag=1
      aqq=poly_fit(alog10(fsmc[0:hiind-1]),alog10(gqqsmc[0:hiind-1]/qq),1)
    ENDIF 
    fend=fq[6:n_elements(fq)-1]
    fendl=alog10(fend)          ;log frequency vector
    gttfit=10^(fendl^3*att[3]+fendl^2*att[2]+fendl*att[1]+att[0])
    if aqq1flag eq 0 then begin
      gqqfit=10^(fendl^2*aqq[2]+fendl*aqq[1]+aqq[0]) 
    endif else begin
      gqqfit=10l^(fendl*aqq[1]+aqq[0]) 
    ENDELSE
    phiflag = 0
    IF hiind GE 6 THEN BEGIN 
      wherephi = where(gqqfit[0:hiind-6] GE gttfit[0:hiind-6],phiflag) ;find freq where gqqfit > gttfit
      IF phiflag GT 0 THEN phistart = wherephi[n_elements(wherephi)-1] ELSE phistart = 0 ;changed from zero on Feb 2003
    ENDIF ELSE phistart = 0
    lofreq=fq[6+phistart]
    IF NOT keyword_set(quiet) THEN print,'    F at phistart: ',lofreq
    qphi1=[[fltarr(1,6+phistart)+1],[transpose(gqqfit[phistart:n_elements(gqqfit)-1]/gttfit[phistart:n_elements(gttfit)-1])]] ;transfer function for phi
    qphi2 = rotate(qphi1,2) ;mirror for other side
    qphi=[[qphi1],[qphi2]]
    tdgq=double(fft(atq*sqrt(qphi),/inverse,/double)) ;degrade time series
    tdgqt=tdgq+tqmn
    power_spectrum,tdgq,hz,/avg,/dtrnd,gxx=gttdg,v=outtdgq,f=fttdg ;plot stats changed 5 to 10
    tdgq = detrend(tdgq)
    tdgq = tdgq - mean(tdgq) ;changed to just taking from v
    gttdgsm=running_avg(cdr(gttdg),11)
    gttdgsmc=gttdgsm[5:n_elements(gttdgsm)-6]
    ttdg=(1/delfq)*total(gttdgsmc[0:14])
    IF NOT keyword_set(noplot) THEN BEGIN
      oldmult = !p.multi
      !p.multi = 0
      subtitle='Spectral Correction'
      IF keyword_set(plottitle) THEN subtitle = plottitle + ' ' + subtitle
      xmin=min(fsmc)
      xmax=max(fend)
      ymin=min([gttsmc/tt,gqqsmc/qq,gttfit,gqqfit])
      ymax=max([gttsmc/tt,gqqsmc/qq,gttfit,gqqfit])
      plot,[xmin,xmax],[ymin,ymax],/nodata,/xlog,/ylog,title=subtitle
      oplot,fsmc,gttsmc/tt
      oplot,fsmc,gttdgsmc/ttdg,linestyle=1
      oplot,fsmc,gqqsmc/qq,color=7600000
      oplot,fend,gttfit,color=7999999
      oplot,fend,gqqfit,color=7999999,linestyle=2
      oplot,[hifreq,hifreq],[1e10,1e-12],linestyle=0
      oplot,[lofreq,lofreq],[1e10,1e-12],linestyle=1
      if badtspecflag eq 1 then begin
         xyouts,0.2,0.15,'T spec < Q spec; no correction done.',/normal
      endif
      if badkneeflag eq 1 then begin
         xyouts,0.2,0.1,'Knee difficulty; no correction done.',/normal
      ENDIF
      !p.multi = oldmult
    ENDIF
    wtq=total(wq*tq)/nq         ;calc correction factor
    wtdgq=total(wq*tdgq)/nq
    q_cfactor=wtq/wtdgq
    qflag=0
    if (q_cfactor gt 2.0) or (q_cfactor lt 1.0) OR (isnan(q_cfactor)) then begin
      qflag=q_cfactor
      q_cfactor=1
    endif
    specs[2:3] = [q_cfactor,qflag]
    if badtspecflag eq 1 then specs[2:3] = [1,-60]
    if badkneeflag eq 1 then specs[2:3] = [1,-50]
  ENDIF

;Find c_cfactor
  IF (clen_ok EQ 1) AND (c_gooddata EQ 1) THEN BEGIN
    nc = 2l^long(alog(lenc)/alog(2))
    dcc = zapbadval(dc[0:nc-1]) & dtcc = zapbadval(dtc[0:nc-1]) & dwcc = zapbadval(dwc[0:nc-1])
    cmn = mean(dcc,/nan) & tcmn = mean(dtcc,/nan) & wcmn = mean(dwcc,/nan)
    c = dcc - cmn & tc = dtcc - tcmn & wc = dwcc - wcmn
    atc=fft(tc,/double) ;Do fourier transform
    halfc=nc/2+1
    fc=dindgen(halfc-1)*fs/(nc-1)
   ;From Lenshow and Raupach phi(f)=exp(-79.7*a*f^2*L*Re^(-1/8)*U^(-2))
   ;where a=tube radius, f=frequency, L=length of tubing,
   ;U=mean flow speed, and Re=2*a*U/nu (Reynolds number) with
   ;nu=molecular viscosity. More simply, phi=exp(el*f^2) where
   ;el=-79.7*a*L*Re^(-1/8)*U^(-2).
    a = co2_const('tubeinnerradius') ;tube inner radius (1/32")
    L = co2_const('tubelength') ;tube length in meters 
    u=L/(lagc+0.0001)
    temp = mean(ndt,/nan)
    nu=(1.18e-4*temp^2+8.844e-2*temp+13.427)*1e-6 ;2nd order fit - Fund. Heat/Mass Trans.
    re=2*a*u/nu                 ;reynolds number
    IF re GT 2300.0 THEN BEGIN
      el=-79.7*a*L*re^(-1/8)*u^(-2)
      IF NOT keyword_set(quiet) THEN print,'    Turbulent tube flow'
    ENDIF ELSE BEGIN
      el=-0.41*a*L*re*nu/(0.16e-4*u^2)
      IF NOT keyword_set(quiet) THEN print,'    Laminar tube flow'
    ENDELSE
    IF NOT keyword_set(quiet) THEN print,'    u = ',u,' temp = ',temp,' nu = ',nu,' re = ',re,' el =',el
    cphi1=exp(el*transpose(fc)^2)
    cphi2 = rotate(cphi1,2)     ;mirror phi 
    cphi = [[cphi1],[cphi2]]  
    tdgc=double(fft(atc*cphi,/inverse,/double))
    tdgct=tdgc+tcmn
    wtc=total(wc*tc)/nc
    wtdgc=total(wc*tdgc)/nc
    c_cfactor=wtc/wtdgc
    cflag=0

    if c_cfactor gt 1.40 or c_cfactor lt 1.0 OR isnan(c_cfactor) then begin
      cflag=c_cfactor
      c_cfactor=1.0
    ENDIF
    specs[0:1] = [c_cfactor,cflag]
  ENDIF 

  IF NOT keyword_set(quiet) THEN print,'    Specs (c,cflag,q,qflag) ',specs
  return,specs
END

PRO save_Spec,specarr,doy,yr
  IF NOT keyword_set(yr) THEN year=co2_const('cur_year') ELSE year = yr
  year2 = year MOD 100
  fhead = co2_const('file_head')
  yrstr = string(Year2,format='(i2.2)')
  specfile = co2_basedir('spec',yr=year)+fhead+yrstr+'.spec.xdf'
  IF file_test(specfile,/read) THEN BEGIN
    specfl = read_xdf(specfile,header=specheader)
  ENDIF ELSE BEGIN
    print,'  Creating SPEC file ',specfile
    specfl = make_array(6,17568,/float,value=nan())
    specheader = ['doy','hhmm','c_cfactor','cflag','q_cfactor','qflag']
    specfl[0,*] =  float(fix(findgen(17568) / 48) + 1)
    specfl[1,*] =  float(fjday_to_hhmm((findgen(17568) MOD 48) / 48.0))
    specfl[2,*] = 1.0
    specfl[3,*] = -40.0
    specfl[4,*] = 1.0
    specfl[5,*] = -40.0
  ENDELSE
  locs = ((doy-1)*48) + hhmm_to_fjday(specarr[0,*]) * 48
  specfl[2:5,locs] = specarr[1:4,*]
  print,'  Writing SPEC file ',specfile
  create_dir,co2_basedir('spec',yr=year)
  write_xdf,specfile,specfl,header=specheader
END

PRO fill_spec,yr,nosave=nosave
;fill missing C data with 45 day time specific mean
;fill missing Q data with 20 day average
;keep only data with flag = 0

  IF NOT keyword_set(yr) THEN year=co2_const('cur_year') ELSE year = yr
  year2 = year MOD 100
  fhead = co2_const('file_head')
  yrstr = string(Year2,format='(i2.2)')
  specfile = co2_basedir('spec',yr=year)+fhead+yrstr+'.spec.xdf'
  IF file_test(specfile,/read) THEN BEGIN
    print,'Reading specfile ',specfile
    specdata = read_xdf(specfile,header=specheader)
  ENDIF ELSE BEGIN
    print,'Cannot find specfile ',specfile
    GOTO,leave
  ENDELSE

;fill c holes
  badspec = where(isnan(specdata[2,*]) OR (specdata[3,*] NE 0.0),nbadspec,complement=goodspec,ncomplement=ngoodspec)
  IF ngoodspec GT 0 THEN wholeyrmean_c = mean(specdata[2,goodspec],/nan) ELSE wholeyrmean_c = 1.0
  IF nbadspec GT 0 THEN BEGIN
    print,'  Filling C spec holes n: ',nbadspec
    specdata[2,badspec] = nan()  ;kill off the bad values
    FOR i = 0,nbadspec-1 DO BEGIN
      badspecloc = badspec[i]
      baddoy = specdata[0,badspecloc]
      badtm = hhmm_to_fjday(specdata[1,badspecloc])
      badtmloc = fix(badtm * 48)
      startdoy = baddoy - 45
      enddoy = baddoy + 45
      IF startdoy LT 1 THEN BEGIN
        enddoy = enddoy - startdoy
        startdoy = 1
      ENDIF
      IF enddoy GT days_in_year(year) THEN BEGIN
        enddoy = days_in_year(year)
        startdoy = startdoy + (enddoy - days_in_year(year))
      ENDIF
      fillloc = ((numgen(startdoy,enddoy) - 1) * 48) + badtmloc
      meanval = nan()
      IF n_elements(specdata[2,fillloc]) GT 7 THEN meanval = mean(specdata[2,fillloc],/nan)
;if can't do time average then create an ensemble interpolated day 45 days
     IF isnan(meanval) THEN BEGIN
        filldata = make_array(48,/float,value=nan())
        FOR x = 0,47 DO BEGIN
          fillloc = ((numgen(startdoy,enddoy) - 1) * 48) + x
          IF n_elements(specdata[2,fillloc]) GT 7 THEN filldata[x] = mean(specdata[2,fillloc],/nan)
        ENDFOR
        filldata = zapbadval(filldata)
        meanval = filldata[badtmloc]
        IF isnan(meanval) THEN meanval = wholeyrmean_c
      ENDIF
      specdata[2,badspecloc] = meanval
      specdata[3,badspecloc] = -70
    ENDFOR
    choles = 1b
  ENDIF ELSE BEGIN
    print,'  No C holes to fill'
    choles = 0b
  ENDELSE

;fill q holes
  badspec = where(isnan(specdata[4,*]) OR (specdata[5,*] NE 0.0),nbadspec,complement=goodspec,ncomplement=ngoodspec)
  IF ngoodspec GT 0 THEN wholeyrmean_q = mean(specdata[4,goodspec],/nan) ELSE wholeyrmean_q = 1.0
  IF nbadspec GT 0 THEN BEGIN
    print,'  Filling Q spec holes n: ',nbadspec
    specdata[4,badspec] = nan() ;kill off the bad values
    startloc = badspec - 960
    endloc = badspec + 960
    badstart = where(startloc LT 0,nbadstart)
    IF nbadstart GT 0 THEN BEGIN
      endloc[badstart] = endloc[badstart] - startloc[badstart]
      startloc[badstart] = 0
    ENDIF
    endpt = days_in_year(year)*48
    badend = where(endloc GE endpt,nbadend)
    IF nbadend GT 0 THEN BEGIN 
      startloc[badend] = startloc[badend] - (endloc[badend] - endpt - 1)
      endloc[badend] = endpt-1
    ENDIF
    meanval = make_array(nbadspec,/float,value=nan())
    FOR i = 0,nbadspec-1 DO BEGIN 
      meanval[i] = mean(specdata[4,startloc[i]:endloc[i]],/nan)
      IF isnan(meanval[i]) THEN meanval[i] = wholeyrmean_q
    ENDFOR
    badmean = where(isnan(meanval),nbadmean)
    IF nbadmean GT 0 THEN BEGIN 
      doys = specdata[0,badspec[badmean]]
      doys = doys[uniq(doys,sort(doys))]
      print,'  Cannot fill holes for all days (DOYS,nvals) ',doys,nbadmean
    ENDIF
    specdata[4,badspec] = meanval
    specdata[5,badspec] = -70.0
    qholes = 1b
  ENDIF ELSE BEGIN
    print,'  No Q holes to fill'
    qholes = 0b
  ENDELSE

  IF ((choles EQ 1) OR (qholes EQ 1)) AND (NOT keyword_set(nosave)) THEN BEGIN
    print,'  Writing cleaned spec file'
    write_xdf,specfile,specdata,header=specheader
  ENDIF

  IF NOT keyword_set(nosave) THEN update_version,'fill_spec',1,year

leave:
  print,'Spec filling complete.'
END