' upper_kuparuk_met.cr1
' by Joel Bailey 7/10/2005
' adapted from Ken Irving's das12.cr1
' constants for base logger-diagnostics precip-te525m rh-pwr-on air-hc2s3 wind-rmyoung air-109 soilsurface-109 barometer-cs106 rh-pwr-off netrad-nrlite snow-sr50a time-variables relay-control-sub camera modem thermistors soilh2o-cs650 pt-aquistar turbidity-obs500 call-met-tables call-hydro-tables call-soil-tables
' 8/18/2015 Joel Bailey - updated Steinhart equations to correct values
' 7/4/2016 Joel Bailey - changed 10m AT to a H2CS3 sensor
' 9/21/2016 Joel Bailey - changed 10m AT to HMP155A sensor
' 3/8/2016 Joel Bailey - changed AT/RH delay from 2 to 4 secs, changed table size of MET table to -1

Const NR_WS_CORR = 0.021286
Const TWOPI = 4 * ATN2(1,0)

Const MUX_PULSE_USEC = 10000

'' Steinhart & Hart constants
'' steinhart-hart-constants for Ground thermistors
Const kA = 0.00128
Const kB = 0.0002368
Const kC = 0.0000000906

'' Steinhart & Hart constants for YSI
'' steinhart-hart-constants for ATB
Const ATB_kA = 0.0014654354
Const ATB_kB = 0.0002386780
Const ATB_kC = 0.0000001000

ConstTable
    Const PAKBUS_INFO = 186
    Const P_MODL = "TE525MM"
    Const P_TIP = 0.1
    Const P_UNITS = "mm"
    Const P__PP = 2
    Const RH_RELAY__CP = 4
    Const RH_1M_MODL = "HP45C"
    Const RH_1M_SN = ""
    Const RH_1M_MULT = 0.1
    Const RH_1M_OFS = 0
    Const AT_1M_MULT = 0.1
    Const AT_1M_OFS = -40
    Const AT_1M__SE = 1
    Const RH_1M__SE = 2
    Const RH_3M_MODL = "HP45C"
    Const RH_3M_SN = ""
    Const RH_3M_MULT = 0.1
    Const RH_3M_OFS = 0
    Const AT_3M_MULT = 0.1
    Const AT_3M_OFS = -40
    Const AT_3M__SE = 3
    Const RH_3M__SE = 4
    Const RH_10M_MODL = "HMP155A"
    Const RH_10M_SN = ""
    Const RH_10M_MULT = 0.1
    Const AT_10M_OFS = -80
    Const AT_10M_MULT = 0.14
    Const RH_10M_OFS = 0
    Const AT_10M__SE = 5
    Const RH_10M__SE = 6
    Const AT_10m__CP = 3
    Const AT_10M__EX = VX3
    Const WS_1M_MODL = "014A MET ONE"
    Const WS_1M_UNITS = "m/s"
    Const WS_1M__PP = 17
    Const WS_3M_MODL = "014A MET ONE"
    Const WS_3M_UNITS = "m/s"
    Const WS_3M__PP = 18
    Const WS_10M_MODL = "RMY 05103"
    Const WS_10M_UNITS = "m/s"
    Const WS_10M_MULT = 0.098
    Const WD_10M_UNITS = "deg"
    Const WD_10M_OFS = 0
    Const WD_10M_MULT = 355
    Const WS_10M__PP = 1
    Const WD_10M__SE = 8
    Const WD_10M__EX = Vx1
    Const ATB_MODL = "GWSCI"
    Const ATB_SN = ""
    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 = 6
    Const NR_MODL = "NR-LITE"
    Const NR_SN = ""
    Const NR_CALIB = 11.9
    Const NR_CALIB_UNITS = "uV/W/m^2"
    Const NR_UNITS = "W/m^2"
    Const NR__DIFF = 5
    Const SNW_MODL = "SR50A"
    Const SNW_SN = ""
    Const SNW_HT_CM = 272.89
    Const SNW_SDI_ADDR = 0
    Const SDI12__CP = 5
    Const MODM_PWR__CP = 9
    Const MODM_MODL = "FGR115"
    Const MODM_SN = ""
    Const MODM_ON_MIN = 0
    Const MODM_OFF_MIN = 10
    Const MODM_PERIOD = 30
    Const MODM_V_LIM = 12.4
    Const MODM_HOLD_V_LIM = 12.6
    Const THERMISTOR_TYPE = ""
    Const THERMISTOR_QTY = 15
    Const MUX_SE_QTY = 1
    Const MUX_THERM__SE = 16
    Const MUX_RESET__CP = 1
    Const MUX_CLOCK__CP = 2
    Const MUX__EX = Vx2
    Const R_REF_KOHMS = 10
    Const HUSKFLUX_QTY = 2
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-hp45c 1m
Public at_1m, rh_1m, dp_1m
Units at_1m = degC
Units dp_1m = degC
Units rh_1m = percent

' air-hp45c 3m
Public at_3m, rh_3m, dp_3m
Units at_3m = degC
Units dp_3m = degC
Units rh_3m = percent

' air-hp35c 10m
Public at_10m, rh_10m, dp_10m
Units at_10m = degC
Units dp_10m = degC
Units rh_10m = percent

' Met_one at 1m
Public ws_1m
Alias ws_1m = windSpeed_1m
Units ws_1m = m/s

' Met_one at 3m
Public ws_3m
Alias ws_3m = windSpeed_3m
Units ws_3m = m/s

' wind-rmyoung at 10m:
Public ws_10m, wd_10m
Alias ws_10m = windSpeed_10m
Units ws_10m = m/s
Alias wd_10m = windDir_10m
Units wd_10m = WD_UNITS

' barometer-cs106:
Public pb
Units pb = mbar

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

' 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(THERMISTOR_QTY), therm_R(THERMISTOR_QTY)
Dim i,therm_X(THERMISTOR_QTY),
Alias therm = st_5cm,st_10cm,st_15cm,st_20cm,st_25cm,st_30cm,st_35cm,st_40cm,st_45cm,st_50cm,st_55cm,st_60cm, radtemp1, radtemp2, atb
Alias therm_R = st_5cm_R,st_10cm_R,st_15cm_R,st_20cm_R,st_25cm_R,st_30cm_R,st_35cm_R,st_40cm_R,st_45cm_R,st_50cm_R,st_55cm_R,st_60cm_R, radtemp1_R, radtemp2_R, atb_R
Units therm_R() = kOhms
Units therm() = degC

' HuskFlux:
Public huskflux(HUSKFLUX_QTY)
Units huskflux() = mV

' 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_1m,fp2,0)
    Average(1,at_3m,fp2,0)
    Average(1,at_10m,fp2,0)
    Average(1,ws_1m,fp2,0)
    Average(1,ws_3m,fp2,0)
    WindVector(1,windSpeed_10m,windDir_10m,fp2,0,0,0,0)
    Average(1,atb,fp2,0)
    Average(1,netRad,fp2,0)
    Average(1,huskflux(),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_1m,fp2)
    Average(1,at_1m,fp2,0)
    Minimum(1,at_1m,fp2,0,0)
    Maximum(1,at_1m,fp2,0,0)
    Sample(1,rh_1m,fp2)
    Minimum(1,rh_1m,fp2,0,0)
    Maximum(1,rh_1m,fp2,0,0)
    Sample(1,dp_1m,fp2)
    Average(1,dp_1m,fp2,0)
    Minimum(1,dp_1m,fp2,0,0)
    Maximum(1,dp_1m,fp2,0,0) 
    Sample(1,at_3m,fp2)
    Average(1,at_3m,fp2,0)
    Minimum(1,at_3m,fp2,0,0)
    Maximum(1,at_3m,fp2,0,0)
    Sample(1,rh_3m,fp2)
    Minimum(1,rh_3m,fp2,0,0)
    Maximum(1,rh_3m,fp2,0,0)
    Sample(1,dp_3m,fp2)
    Average(1,dp_3m,fp2,0)
    Minimum(1,dp_3m,fp2,0,0)
    Maximum(1,dp_3m,fp2,0,0)
    Sample(1,at_10m,fp2)
    Average(1,at_10m,fp2,0)
    Minimum(1,at_10m,fp2,0,0)
    Maximum(1,at_10m,fp2,0,0)
    Sample(1,rh_10m,fp2)
    Minimum(1,rh_10m,fp2,0,0)
    Maximum(1,rh_10m,fp2,0,0)
    Sample(1,dp_10m,fp2)
    Average(1,dp_10m,fp2,0)
    Minimum(1,dp_10m,fp2,0,0)
    Maximum(1,dp_10m,fp2,0,0)
    Sample(1,windSpeed_1m,fp2)
    Sample(1,windSpeed_3m,fp2)
    Sample(1,windSpeed_10m,fp2)
    Sample(1,windDir_10m,fp2)
    WindVector(1,windSpeed_10m,windDir_10m,fp2,0,0,0,0)
    Maximum(1,windSpeed_1m,fp2,0,0)
    Maximum(1,windSpeed_3m,fp2,0,0)
    Maximum(1,windSpeed_10m,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,huskflux(),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

' HrlySoil
DataTable(HrlySoil,1,-1)
    DataInterval(0,60,Min,0)
    Sample(1,progSig,uint2)
    Average(THERMISTOR_QTY,therm(),fp2,0)
EndTable

' Met: base precip-te525m air-hc2s3 wind-rmyoung air-109 barometer-cs106 netrad-nrlite snow-sr50a 
DataTable(Met,1,-1)
    DataInterval(0,60,Min,0)
    Sample(1,progSig,uint2)
    Totalize(1,precip,fp2,0)
    Average(1,at_1m,fp2,0)
    Sample(1,rh_1m,fp2)
    Average(1,at_3m,fp2,0)
    Sample(1,rh_3m,fp2)
    Average(1,at_10m,fp2,0)
    Sample(1,rh_10m,fp2)
    Average(1,windSpeed_1m,fp2,0)
    Average(1,windSpeed_3m,fp2,0)
    WindVector(1,windSpeed_10m,windDir_10m,fp2,0,0,0,0)
    Average(1,atb,fp2,0)
    Average(1,pb,fp2,0)
    Average(1,netRad,fp2,0)
    Average(1,huskflux(),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_3m,fp2,0)
    Average(1,snowdepth,fp2,0)
    Sample(1,sr50_height,fp2)
    Average(1,sr50_quality,fp2,0)
EndTable

' Therm-R: base thermistors 
DataTable(thermRes,1,-1)
    DataInterval(0,60,Min,0)
    Sample(1,progSig,uint2)
    Average(THERMISTOR_QTY,therm_R(),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 = 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 j, lnR
    For j=1 To THERMISTOR_QTY
        If j < 15 Then
          '' calculate thermistor resistance from ratio and ref R:
          therm_R(j) = R_REF_KOHMS * therm_X(j) / (1-therm_X(j))
          '' calculate natural log of thermistor resistance:
          lnR = LN(therm_R(j)*1000)
          '' use Steinhart & Hart equation to derive temperature:
          therm(j) = (1/(kA + kB*lnR + kC*lnR^3)) - 273.15
        Else ' calculates ATB thermistor based on different Steinhart constants
          '' calculate thermistor resistance from ratio and ref R:
          therm_R(j) = R_REF_KOHMS * therm_X(j) / (1-therm_X(j))
          '' calculate natural log of thermistor resistance:
          lnR = LN(therm_R(j)*1000)
          '' use Steinhart & Hart equation to derive temperature:
          therm(j) = (1/(ATB_kA + ATB_kB*lnR + ATB_kC*lnR^3)) - 273.15
        EndIf  
    Next j
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,4000,mSec)
    
        ' air-hp45c 1m:
        VoltSe(at_1m,1,mV2500,AT_1M__SE,1,0,250,AT_1M_MULT,AT_1M_OFS)
        VoltSe(rh_1m,1,mV2500,RH_1M__SE,1,0,250,RH_1M_MULT,RH_1M_OFS)
        DewPoint(dp_1m,at_1m,rh_1m)
        
       ' air-hp45c 3m:
        VoltSe(at_3m,1,mV2500,AT_3M__SE,1,0,250,AT_3M_MULT,AT_3M_OFS)
        VoltSe(rh_3m,1,mV2500,RH_3M__SE,1,0,250,RH_3M_MULT,RH_3M_OFS)
        DewPoint(dp_3m,at_3m,rh_3m)
        
       ' air-mp155c 10m:
        VoltSe(at_10m,1,mV2500,AT_10M__SE,1,0,250,AT_10M_MULT,AT_10M_OFS)
        VoltSe(rh_10m,1,mV2500,RH_10M__SE,1,0,250,RH_10M_MULT,RH_10M_OFS)
        DewPoint(dp_10m,at_10m,rh_10m)
    
		    '014A Wind Speed Sensor measurement 'WS_1m'
		    PulseCount(ws_1m,1,WS_1M__PP,2,1,0.8,0.447)
		    If ws_1m<0.457 Then ws_1m=0
		
        '014A Wind Speed Sensor measurement 'WS_3m'
		    PulseCount(ws_3m,1,WS_3M__PP,2,1,0.8,0.447)
		    If ws_3m<0.457 Then ws_3m=0

        ' wind-rmyoung 10m:
        PulseCount(ws_10m,1,WS_10M__PP,1,1,WS_10M_MULT,0)
        BrHalf(wd_10m,1,mV2500,WD_10M__SE,WD_10M__EX,1,2500,1,0,250,WD_10M_MULT,WD_10M_OFS)
    
        ' barometer-cs106:
        PortSet(PB__CP,1)
        Delay(0,150,mSec)
        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_10m > 5 Then netRad = netRad_raw * (1 + NR_WS_CORR * (windSpeed_10m-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_3m+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)
        ' thermistors:
        PortSet(MUX_RESET__CP,1)
        i = 1
        SubScan(0,Sec,THERMISTOR_QTY)
            PulsePort(MUX_CLOCK__CP,MUX_PULSE_USEC)
            BrHalf (therm_X(i),1,mV2500,MUX_THERM__SE,MUX__EX,1,2500,True,0,250,1.0,0)
            i = i + 1
        NextSubScan
        PortSet(MUX_RESET__CP,0)
        calc_thermistors()
    
        ' call-soil-tables:
        calltable HrlySoil
        CallTable ThermRes
    
    NextScan

endProg