' imnavait_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/26/2015 updated ATB constants, changed steinhart equations for soil thermistors
' 11/18/2015 correctd the stein-hart equations for the soil temp sensors, chaged resistance values to record Kohms instead of ohms

Const TWOPI = 4 * ATN2(1,0)

Const MUX_PULSE_USEC = 10000

'' Steinhart & Hart constants for soil temp sensors
Const kA = 0.0014654354
Const kB = 0.0002386780
Const kC = 0.0000001000

'' 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 = 185
    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 = "HC2S3"
    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 = "HC2S3"
    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 = "HC2S3"
    Const RH_10M_SN = ""
    Const RH_10M_MULT = 0.1
    Const RH_10M_OFS = 0
    Const AT_10M_MULT = 0.1
    Const AT_10M_OFS = -40
    Const AT_10M__SE = 5
    Const RH_10M__SE = 6
    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 = "ATB_GW_SCI"
    Const ATB_SN = ""
    Const ATB_MULT = 1
    Const ATB_OFS = 0
    Const ATB__SE = 7
    Const ATB__EX = Vx3
    Const ATB_REF_OHMS = 1000
    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 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 = 5
    Const SNW_MODL = "SR50A"
    Const SNW_SN = "4415"
    Const SNW_HT_CM = 176
    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 = 14
    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 = 4
    Const HUSKFLUX_QTY = 2
    Const MUX2_RESET__CP = 3
    Const MUX2_CLOCK__CP = 6
    Const MUX2__DIFF = 6
EndConstTable

' base:
Public progSig As Long

' logger-diagnostics:
public battery, SolarPnlV, 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-hp45c 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

' atb GWS
Public atb, atb_X, atb_R
Units atb = degC
Units atb_X = volts
Units atb_R = ohms

' 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 = nc_1,nc_2,st_plus_5cm_snow,st_sfc,st_5cm,st_10cm,st_20cm,st_30cm,st_40cm,st_50cm,st_60cm,st_80cm, st_100cm, st_120cm
Alias therm_R = nc_1_R,nc_2_R,st_plus_5cm_snow_R,st_sfc_R,st_5cm_R,st_10cm_R,st_20cm_R,st_30cm_R,st_40cm_R,st_50cm_R,st_60cm_R,st_80cm_R, st_100cm_R, st_120cm_R
Units therm_R() = KOhms
Units therm() = degC

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

' 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: base soilsurface-109 thermistors soilh2o-cs650 
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,8)
    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(therm_R,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
        '' 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
    Next j
EndSub

' thermistors:calc_atb()
Sub calc_atb()
    Dim lnR
        '' calculate thermistor resistance from ratio and ref R:
        atb_R = ATB_REF_OHMS * atb_X / (1-atb_X)
        '' calculate natural log of thermistor resistance:
        lnR = LN(atb_R)
        '' use Steinhart & Hart equation to derive temperature:
        atb = (1/(ATB_kA + ATB_kB*lnR + ATB_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-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-hp45c 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)
    
        ' atb GWS
        BrHalf (atb_X,1,mV2500,ATB__SE,ATB__EX,1,2500,True,0,250,ATB_MULT,ATB_OFS)
        calc_atb()
    
        ' barometer-cs106:
        VoltSe(pb,1,mV2500,PB__SE,1,0,250,PB_MULT,PB_OFS)
    
        ' 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()
        
        'Huskflux radiometers
        PortSet(MUX2_RESET__CP,1)
        j=1
        SubScan(0,Sec,HUSKFLUX_QTY)
          ' Husk Flux on AM416:
          PulsePort(MUX2_CLOCK__CP,MUX_PULSE_USEC)
          VoltDiff (huskflux(i),1,mV2500,MUX2__DIFF,True,0,250,1.0,0)
          j = j + 1
        NextSubScan
        PortSet(MUX2_RESET__CP,0)
    
        ' call-soil-tables:
        calltable HrlySoil
        calltable Therm_R
    
    nextScan

endProg