! 31 Dec 05

  SysCmd 'echo ProcEppleys >con'

  ! Upward facing pyrgeometer (read as volts)
  Td = {lUPTd0}+lUPdme*({lUPTd1}+lUPdme*({lUPTd2}+lUPdme*({lUPTd3}+lUP*{lUPTd4})))
  Td = Td+273.2
  Ts = {lUPTs0}+lUPbdy*({lUPTs1}+lUPbdy*({lUPTs2}+lUPbdy*({lUPTs3}+lUPbdy*{lUPTs4})))
  Ts = Ts+273.2
  Raw = {lUPTh1}*lUP+{lUPTh0}
  lUP = ((Raw-{lUPThO})/{lUPThk})-4.*SB*Td^4+5.*SB*Ts^4
  Clear lUPdme lUPbdy

  ! Downward facing pyrgeometer (read as Volts)
  Td = {lDNTd0}+lDNdme*({lDNTd1}+lDNdme*({lDNTd2}+lDNdme*({lDNTd3}+lDNdme*{lDNTd4})))
  Td = Td+273.2
  Ts = {lDNTs0}+lDNbdy*({lDNTs1}+lDNbdy*({lDNTs2}+lDNbdy*({lDNTs3}+lDNbdy*{lDNTs4})))
  Ts = Ts+273.2
  Raw = {lDNTh1}*lDN+{lDNTh0}
  lDN = ((Raw-{lDNThO})/{lDNThk})-4.*SB*Td^4+5.*SB*Ts^4
  Clear lDNfme lDNbdy

  ! Eppley downward facing pyranometer.
  sDN = sDN * {sDNCal}

  ! Get the raw incoming shortwave radiation.
  sUP = sUP * {sUPCal}
  sUPraw = sUP

  ps = ps_RB
  taRF = tbrec_RB

  $IF {sUPCor} $EQ Y $THEN

    ! Correct the incoming shortwave for the radiometer tilt with respect to
    ! the aircraft and the aircraft's attitude with respect to the sun.
    ! The method follows Katsaros & DeVault (1986) but uses the expression for
    ! the incidence angle given in Saunders et al (1991).
    ! First, we get the solar azimuth and zenith angles from the date, time,
    ! latitude and longitude.
   
    Almanac az zd zdr ref Rlat Rlon UTCtime UTCdate ps taRF
    Clear zdr ref
   
    ! Filter the aircraft attitude data to better match the response time
    ! of the Eppley PSP (about 1 second) and the AHRS.  We use a simple exponential
    ! filter (non-symmetric) and choose the filter response to minimise the
    ! phase difference between the shortwave and attitude data.
   
    Hdg = Rthdg
    Pch = Rpch
    Rll = Rrll
   
    DirnHelix Hdg Hdg 270

    FILTER
     FilterType NOAA
     DataInt    dt
     Pass       Low 1.0
     Select     Hdg Pch Rll
    END

    CreateLs  sdmy  Hdg  1.0  1.0		! Create a series of "1.0"
    u,v(gs,hdg) HdgX HdgY sdmy  Hdg		! Split ahdg into components
    gs,hdg(u,v) sdmy  Hdg  HdgX HdgY		! Filtered ahdg from components
    Clear sdmy
   
    ! Add the pitch and roll offsets of the radiomater mounting.
    sUPPch  = Pch + sUPPOff
    sUPRll  = Rll + sUPROff
    Clear Pch Rll

    ! Get the terms for the attitude correction.
    spch    = sin(sUPPch)	! Sine of the radiometer pitch angle
    cpch    = cos(sUPPch)	! Cosine of the radiometer pitch angle
    srll    = sin(sUPRll)	! Sine of the radiometer roll angle
    crll    = cos(sUPRll)	! Cosine of the radiometer roll angle
    shdg_az = sin(Hdg-az)	! Sine of heading minus solar azimuth
    chdg_az = cos(Hdg-az)	! Cosine of heading minus solar azimuth
    szd     = sin(zd)		! Sine of solar zenith
    czd     = cos(zd)           ! Cosine of solar zenith
   
    ! We are ready to calculate the incidence angle, the pyranometer calibration
    ! and the attitude correction.
    ! First, we calculate the cosine of the incidence angle, the angle
    ! between the direct beam and the normal to the thermopile.  The following
    ! expression is from Saunders et al (1991).
   
    cdlt = -crll*spch*szd*chdg_az-srll*szd*shdg_az+crll*cpch*czd
   
    ! Now get the correction factor, R, but we trap times when the solar zenith
    ! is greater than 80 deg ie czd <= 0.17 and set R to 1.0 to disable the
    ! attitude correction.
   
    Div0 R cdlt czd 1.0 0.17
   
    ! Next, we limit the incidence angle to values less than 80 deg ie
    ! cdlt >= 0.17 and interpolate R when the incidence angle is greater
    ! than this.  This is done to avoid spikes in the corrected incoming
    ! shortwave at high incidence angles.
   
    PICK
      Indicator cdlt
      Between   0.17 1.0
      Inpol     linear
      Select    R
    END
   
    ! We now need to decide when to apply the correction based on whether the
    ! incoming shortwave is dominated by direct or diffuse beam radiation.
    ! First, we remove low sUP values due to non-zero aircraft pitch and
    ! roll, this is done to make sure that these are not mistaken as cloud
    ! shadows.  We remove the low values by rejecting data when the incidence
    ! angle exceeds +/- 10 deg of the solar zenith and replacing the rejected
    ! data with interpolated values.
   
    DUpr = cos(zd-5.) - cdlt
    DLwr = cdlt - cos(zd+5.)
    CondSampInd LTUpr DUpr dt 0 999999 1
    CondSampInd GTLwr DLwr dt 0 999999 1
    PInd = LTUpr + GTLwr
    PICK
      Indicator PInd
      Between   1.5 2.5
      Inpol     linear
      Select    sUP
    END
    Clear LTUpr GTLwr DUpr DLwr PInd
   
    ! Now we calculate the threshold value used to decide when diffuse dominates
    ! over direct.  Note that Saunders et al (1991) use Scrit = 920.*czd^1.28
   
    Scrit = 800.*(czd^1.28)	! Threshold for diffuse/direct
   
    ! Next, we construct two series, one called GEScrit which is 1.0 when sUP
    ! is greater than Scrit and the other called LTScrit which is 1.0 when
    ! sUP is less than Scrit.  Note that periods of sUP < Scrit shorter than
    ! 2.5 seconds (~100 m) are ignored.
   
    SD = sUP - Scrit
    CondSampInd GEScrit SD dt 0 999999 2.5	! 1.0 when sUP >= Scrit
    InvertCSI   LTScrit GEScrit		! 1.0 when sUP < Scrit
   
    ! Now apply the correction as required.
   
    R = LTScrit + GEScrit*({Eps}*R+(1.-{Eps}))
    sUP = sUPraw/R
    Clear spch cpch srll crll shdg_az chdg_az szd SD GEScrit LTScrit
  $ENDIF
   
  $IF Y $IN {PltAnaEppleys} $THEN
    abscissa = GPSsec
    PLOT
      $LOOP  ser1    sUP   sDN   lUP   lDN
      $&     ser2    none  none none  none  
      $&     ser3    none  none none  none  
      $&     ser4    none  none none  none  
      $&     ser5    none  none none  none  
      $&     pltdev  N     N    N     N     
      $INDEX k
        $IF {k} $EQ 1 $THEN
          Size {Orient}
          Logo off
          NewRow   1
          Header * '{GRP} {EXPT} {AC} {run} {Date} GPSsec: {Gstart}-{Gend}'
          Header * '   '
          Header * '   '
        $ENDIF
        $INCLUDE {CIncPath}procs\newplot5mnNH.inc
      $ENDLOOP
    END
  $ENDIF