' put_met.cri   
' adapted from kcirving program das2-12sep2012.cri, modified by Joel Bailey  7/8/2015
' constants for base logger-diagnostics precip-te525m rh-pwr-on air-hp45c wind-rmyoung air-109 barometer-cs106 rh-pwr-off netrad-nrlite snow-sr50a time-variables relay-control-sub thermistors call-met-tables
' 7/29/2015 - Joel Bailey - Changed on/off time of modem to once an hour instead of every 30 mins
' 7/31/2015 - Joel Bailey - Added 2 second delay on CP8 for PB measurement, changed Steinhart Constants for ATB, added RadRaw readings
' 5/27/2016 - Joel Bailey - added individual constants to Eppley Radiometers, adjusted for correct orientation of UP and Down per Rob Gieck
' 7/19/2016 - Joel Bailey - changed Raw data
' 7/21/2016 - Joel Bailey - added temperature measurements on all 4 channels of Multiplexor, in case where LW and SW are not plugged into same channels every time, added more to Radiation raw data table
' 7/26/2016 - Joel Bailey - swtiched Raw data resistance to K ohms
Const NR_WS_CORR = 0.021286

'' Steinhart & Hart constants
Const kA = 0.0014654354
Const kB = 0.0002386780
Const kC = 0.0000001000
Const pir_kA = 0.0010295
Const pir_kB = 0.0002391
Const pir_kC = 0.0000001568

ConstTable
    Const PAKBUS_INFO = 180
    Const P_MODL = "TE525MM"
    Const P_TIP = 0.1
    Const P_UNITS = "mm"
    Const P__PP = 2
    Const RH_RELAY__CP = 4
    Const RH_MODL = "HP45C"
    Const RH_SN = ""
    Const RH_MULT = 0.1
    Const RH_OFS = 0
    Const AT_MULT = 0.1
    Const AT_OFS = -40
    Const AT__SE = 1
    Const RH__SE = 2
    Const WS_MODL = "RMY 05103"
    Const WS_UNITS = "m/s"
    Const WS_MULT = 0.098
    Const WD_UNITS = "deg"
    Const WD_OFS = 0
    Const WD_MULT = 355
    Const WS__PP = 1
    Const WD__SE = 8
    Const WD__EX = Vx1
    Const ATB_MODL = "ATB_GW_SCI"
    Const ATB_SN = ""
    Const ATB_MULT = 1
    Const ATB_OFS = 0
    Const ATB__SE = 5
    Const ATB__EX = Vx3
    Const PB_MODL = "CS106"
    Const PB_SN = " "
    Const PB_UNITS = "mBar"
    Const PB_MULT = 0.24
    Const PB_OFS = 600
    Const PB__SE = 15
    Const PB__CP = 8
    Const NR_MODL = "NR-LITE"
    Const NR_SN = "114128"
    Const NR_CALIB = 11.9
    Const NR_CALIB_UNITS = "uV/W/m^2"
    Const NR_UNITS = "W/m^2"
    Const NR__DIFF = 2
    Const PIR_MODL_LW_UP = "Eppley PIR"
    Const PIR_SN_LW_UP = "25982F3"
    Const PIR_MULT_LW_UP = 227.273
    Const PIR_MODL_LW_DOWN = "Eppley PIR"
    Const PIR_SN_LW_DOWN = "27597F3"
    Const PIR_MULT_LW_DOWN = 275.482
    Const PSP_MODL_SW_UP = "Eppley PSP"
    Const PSP_SN_SW_UP = "16750F3"
    Const PSP_MULT_SW_UP = 106.609
    Const PSP_MODL_SW_DOWN = "Eppley PSP"
    Const PSP_SN_SW_DOWN = "30156F3"
    Const PSP_MULT_SW_DOWN =124.224
    Const RAD__DIFF = 5
    Const RAD_TEMP__SE = 11
    Const SNW_MODL = "SR50A"
    Const SNW_SN = "4415"
    Const SNW_HT_CM = 152
    Const SNW_SDI_ADDR = 0
    Const SDI12__CP = 5
    Const MODM_PWR__CP = 9
    Const MODM_MODL = "9522B"
    Const MODM_SN = ""
    Const MODM_ON_MIN = 0
    Const MODM_OFF_MIN = 10
    Const MODM_PERIOD = 60
    Const MODM_V_LIM = 12.4
    Const MODM_HOLD_V_LIM = 12.6
    Const THERMISTOR_TYPE = "PS302J2"
    Const THERMISTOR_QTY = 1
    Const R_REF_OHMS = 1000
    Const MUX_RESET__CP = 1
    Const MUX_CLOCK__CP = 2
    Const NO_OF_EPPLEY = 4
EndConstTable

' base:
Public progSig As Long

' logger-diagnostics:
Public battery, PanelT
Units PanelT = degC
Units battery = Volts

' precip-te525m:
Public precip
Units precip = P_UNITS

' air-hc2s3:
Public at, rh, dp
Units at = degC
Units dp = degC
Units rh = percent

' wind-rmyoung:
Public ws, wd
Alias ws = windSpeed
Units ws = WS_UNITS
Alias wd = windDir
Units wd = WD_UNITS

' air-109:
Public atb
Units atb = degC

' barometer-cs106:
Public pb
Units pb = mbar

' netrad-nrlite:
Public nr(2)
Alias nr = netRad_raw,netRad
Units nr() = NR_UNITS

' eppley-pir
Public rad_mv(4), rad_therm_X(4), rad_therm_R(4), rad_temp(4), rad(4)
Alias rad = rad_lw_up, rad_lw_down, rad_sw_up, rad_sw_down
Units rad_mv() = mV
Units rad_therm_X() = mV
Units rad_therm_R() = Kohms
Units rad_temp() = degK
Units rad() = w/m2

' snow-sr50a:
Public sr50_raw(2), sr50_down, snowdepth, sr50_height
Alias sr50_raw(2) = sr50_quality

' time-variables:
Dim Time(9)
Alias Time(4) = thisHour
Alias Time(5) = thisMinute
Alias Time(9) = thisDOY

' modem:
Public modemOnMin, modemOffMin, modemPeriod, modemVlim, modemHldVlim
Public modemHold As Boolean, modemForce As Boolean

' thermistors:
Public therm_R
Dim therm_X
Units therm_R = Ohms

' Daily: base precip-te525m air-hc2s3 wind-rmyoung air-109 netrad-nrlite snow-sr50a 
DataTable(Daily,1,-1)
    DataInterval(0,1440,Min,0)
    Sample(1,progSig,uint2)
    Totalize(1,precip,fp2,0)
    Average(1,at,fp2,0)
    WindVector(1,windSpeed,windDir,fp2,0,0,0,0)
    Average(1,atb,fp2,0)
    Average(1,netRad,fp2,0)
    Average(1,rad_lw_up,fp2,0)
    Average(1,rad_lw_down,fp2,0)
    Average(1,rad_sw_up,fp2,0)
    Average(1,rad_sw_down,fp2,0)
    Average(1,snowdepth,fp2,0)
EndTable

' HrlyAtms: base precip-te525m air-hc2s3 wind-rmyoung air-109 barometer-cs106 netrad-nrlite snow-sr50a 
DataTable(HrlyAtms,1,-1)
    DataInterval(0,60,Min,0)
    Sample(1,progSig,uint2)
    Totalize(1,precip,fp2,0)
    Sample(1,at,fp2)
    Average(1,at,fp2,0)
    Minimum(1,at,fp2,0,0)
    Maximum(1,at,fp2,0,0)
    Sample(1,rh,fp2)
    Minimum(1,rh,fp2,0,0)
    Maximum(1,rh,fp2,0,0)
    Sample(1,dp,fp2)
    Average(1,dp,fp2,0)
    Minimum(1,dp,fp2,0,0)
    Maximum(1,dp,fp2,0,0)
    Sample(1,windSpeed,fp2)
    Sample(1,windDir,fp2)
    WindVector(1,windSpeed,windDir,fp2,0,0,0,0)
    Maximum(1,windSpeed,fp2,0,0)
    Sample(1,atb,fp2)
    Average(1,atb,fp2,0)
    Minimum(1,atb,fp2,0,0)
    Maximum(1,atb,fp2,0,0)
    Sample(1,pb,fp2)
    Average(1,pb,fp2,0)
    Average(1,netRad,fp2,0)
    Average(1,rad_lw_up,fp2,0)
    Average(1,rad_lw_down,fp2,0)
    Average(1,rad_sw_up,fp2,0)
    Average(1,rad_sw_down,fp2,0)
    Average(1,snowdepth,fp2,0)
EndTable

' HrlyDiag: base logger-diagnostics 
DataTable(HrlyDiag,1,-1)
    DataInterval(0,60,Min,0)
    Sample(1,progSig,uint2)
    Average(1,battery,fp2,0)
    Average(1,PanelT,fp2,0)
EndTable

' Met: base precip-te525m air-hc2s3 wind-rmyoung air-109 barometer-cs106 netrad-nrlite snow-sr50a 
DataTable(Met,1,8)
    DataInterval(0,60,Min,0)
    Sample(1,progSig,uint2)
    Totalize(1,precip,fp2,0)
    Average(1,at,fp2,0)
    Sample(1,rh,fp2)
    WindVector(1,windSpeed,windDir,fp2,0,0,0,0)
    Average(1,atb,fp2,0)
    Average(1,pb,fp2,0)
    Average(1,netRad,fp2,0)
    Average(1,rad_lw_up,fp2,0)
    Average(1,rad_lw_down,fp2,0)
    Average(1,rad_sw_up,fp2,0)
    Average(1,rad_sw_down,fp2,0)
    Average(1,snowdepth,fp2,0)
EndTable

' SR50: base air-hc2s3 snow-sr50a 
DataTable(SR50,1,-1)
    DataInterval(0,60,Min,0)
    Sample(1,progSig,uint2)
    Average(1,at,fp2,0)
    Average(1,snowdepth,fp2,0)
    Sample(1,sr50_height,fp2)
    Average(1,sr50_quality,fp2,0)
EndTable


' Epply Radiometer RAW: 
DataTable(RadRaw,1,-1)
    DataInterval(0,60,Min,0)
    Sample(1,progSig,uint2)
    Average(1,rad_mv(1),fp2,0)
    Average(1,rad_mv(2),fp2,0)
    Average(1,rad_mv(3),fp2,0)
    Average(1,rad_mv(4),fp2,0)
    Average(1,rad_therm_X(1),fp2,0)
    Average(1,rad_therm_X(2),fp2,0)
    Average(1,rad_therm_X(3),fp2,0)
    Average(1,rad_therm_X(4),fp2,0)
    Average(4,rad_therm_R(),fp2,0)
    Average(4,rad_temp(),fp2,0)
EndTable

' relay-control-sub:relayControl()
Sub relayControl(relayPort,onMin,offMin,period,Vlimit,holdVlimit,hold,force)
    ' uses globals thisMinute, battery
    If force Then Call setRelayState(relayPort,1)
    If thisMinute MOD period = onMin AND battery >= Vlimit Then Call setRelayState(relayPort,1)
    If thisMinute MOD period = onMin AND battery < holdVlimit Then hold = 0
    If thisMinute MOD period = offMin AND NOT hold Then Call setRelayState(relayPort,0)
EndSub

' relay-control-sub:setPortState()
Sub setRelayState(relay,state As Boolean)
    ' portSet() does not allow variable for port, so brute force it:
    If relay = 9 Then
        PortSet(9,state)
    Else If relay = 8 Then
        PortSet(8,state)
    Else If relay = 7 Then
        PortSet(7,state)
    Else If relay = 6 Then
        PortSet(6,state)
    Else If relay = 5 Then
        PortSet(5,state)
    Else If relay = 4 Then
        PortSet(4,state)
    Else If relay = 3 Then
        PortSet(3,state)
    Else If relay = 2 Then
        PortSet(2,state)
    Else If relay = 1 Then
        PortSet(1,state)
    EndIf
EndSub


' thermistors:calc_thermistors()
Sub calc_thermistors()
    Dim lnR
        '' calculate thermistor resistance from ratio and ref R:
        therm_R = R_REF_OHMS * therm_X / (1-therm_X)
        '' calculate natural log of thermistor resistance:
        lnR = LN(therm_R)
        '' use Steinhart & Hart equation to derive temperature:
        atb = (1/(kA + kB*lnR + kC*lnR^3)) - 273.15
EndSub
BeginProg

    ' for base:
    progSig = status.ProgSignature
    
    ' for snow-sr50a:
    sr50_height = SNW_HT_CM
    
    ' for modem:
    modemOnMin = (MODM_ON_MIN + MODM_PERIOD) MOD MODM_PERIOD
    modemOffMin = (MODM_OFF_MIN + MODM_PERIOD) MOD MODM_PERIOD
    modemPeriod = MODM_PERIOD
    modemVlim = MODM_V_LIM
    modemHldVlim = MODM_HOLD_V_LIM
    
    Scan(60,sec,0,0)
        ' logger-diagnostics:
        PanelTemp(PanelT,250)
        battery(battery)
     
        ' precip-te525m:
        PulseCount(precip,1,P__PP,2,0,P_TIP,0)
    
        ' rh-pwr-on:
        PortSet(RH_RELAY__CP,1)
        Delay(0,2000,mSec)
    
        ' air-hc2s3:
        VoltSe(at,1,mV2500,AT__SE,1,0,250,AT_MULT,AT_OFS)
        VoltSe(rh,1,mV2500,RH__SE,1,0,250,RH_MULT,RH_OFS)
        DewPoint(dp,at,rh)
    
        ' wind-rmyoung:
        PulseCount(ws,1,WS__PP,1,1,WS_MULT,0)
        BrHalf(wd,1,mV2500,WD__SE,WD__EX,1,2500,1,0,250,WD_MULT,WD_OFS)
    
        ' air-109:
        'Therm109(atb,1,ATB__SE,ATB__EX,0,250,ATB_MULT,ATB_OFS)
        BrHalf (therm_X,1,mV2500,ATB__SE,ATB__EX,1,2500,True,0,250,ATB_MULT,ATB_OFS)
        
		    'run subroutine to calculate degrees for each thermistor
		    calc_thermistors()
        
        ' barometer-cs106: allow 2 secs for warm up
        PortSet(PB__CP,1)
        Delay (0,2,Sec)
        VoltSe(pb,1,mV2500,PB__SE,1,0,250,PB_MULT,PB_OFS)
        PortSet(PB__CP,0)
    
        ' rh-pwr-off:
        PortSet(RH_RELAY__CP,0)
    
        ' netrad-nrlite:
        VoltDiff(nr,1,mv25,NR__DIFF,True,0,250,1000/NR_CALIB,0)
        netRad = netRad_raw
        If windSpeed > 5 Then netRad = netRad_raw * (1 + NR_WS_CORR * (windSpeed-5))
        
        ' snow-sr50a:
        '' read SR50A snow depth sensor:
        ''    read distance in m (positive) and quality number into sr50_raw()
        ''    make distance down negative, convert to cm
        ''    apply correction for speed of sound as f(temperature, 0C)
        ''    calculate snow depth as sum of sensor height and distance down
        '' 
        SDI12Recorder(sr50_raw, SDI12__CP, SNW_SDI_ADDR, "M1!", 1, 0)
        sr50_raw = sr50_raw * -100.0
        sr50_down = sr50_raw * SQR((at+273.15) / 273.15)
        snowdepth = sr50_height + sr50_down
    
        ' time-variables:
        RealTime(Time())
    
        ' modem:
        relayControl(MODM_PWR__CP,modemOnMin,modemOffMin,modemPeriod,modemVlim,modemHldVlim,modemHold,modemForce)
    
        ' call-met-tables:
        CallTable Daily
        CallTable HrlyAtms
        CallTable HrlyDiag
        CallTable SR50
        CallTable Met    
    NextScan
    
SlowSequence
	Scan(600,sec,0,0)	  
	'Turn AM16/32 Multiplexer On
		PortSet(MUX_RESET__CP,1)
		Dim i
		i=1
		SubScan(0,uSec,NO_OF_EPPLEY)
			'Switch to next AM16/32 Multiplexer channel
			PulsePort(MUX_CLOCK__CP,10000)
			' read 2 PIR with thermistors, read 2 PSP
      VoltDiff(rad_mv(i),1,mv25,RAD__DIFF,True,0,250,1,0)
      BrHalf (rad_therm_X(i),1,mV2500,RAD_TEMP__SE,Vx2,1,2500,True ,0,250,1.0,0)
      Dim lnR
      '' calculate PIR thermistor resistance from ratio and ref R:
      rad_therm_R(i) = 1 * (rad_therm_X(i) / (1-rad_therm_X(i)))
      '' calculate natural log of thermistor resistance:
      lnR = LN(rad_therm_R(i)*1000)
      '' use Steinhart & Hart equation to derive temperature:
      rad_temp(i) = (1/(pir_kA + pir_kB*lnR + pir_kC*lnR^3))
			i=i+1
		NextSubScan
		rad(1) = (rad_mv(1)*PIR_MULT_LW_UP) + (5.6697E-8 * rad_temp(1)^4)
		rad(2) = (rad_mv(2)*PIR_MULT_LW_DOWN) + (5.6697E-8 * rad_temp(2)^4)
		rad(3) = rad_mv(3)* PSP_MULT_SW_UP
		rad(4) = rad_mv(4)* PSP_MULT_SW_DOWN
		'Turn AM16/32 Multiplexer Off
		PortSet(MUX_RESET__CP,0)
					
    'call rad tables
     CallTable Daily
     CallTable HrlyAtms
     CallTable Met 
     CallTable RadRaw
  NextScan
EndProg