' WERC FB met 10m cr1000 datalogger program -- jb,ki 15sep2010
' 9-15-2016 - Changed AT_10m constants to reflect HPM155 sensor coefficients
' 9-20-2016 - Changed Net_mv to Net_rad, added the constant in the program so records in wpm2 and not mV, added atb 107 at 1m
' 7/8/2017 eky changed themp 107 for atb from SE9 to SE8
StationName(FB_met)
Const PAKBUS_ID = 2

'' DECLARE CHANGEABLE CONSTANTS:
'' NOTE: only first 7 characters are visible on keypad
'' WARNING: changes WILL restart program, MAY lose data!!

ConstTable

    '' Sonic Ranging snow depth sensor
    const SR50_MDL = "SR50A"
    const SR50_SN = "2049"
    const SR50_SDI_ADDR = 0   
    ''   sensor to ground in cm:
    const SR50_HT_CM = 193
    const SR50_AT_SELECT = 2 '' choose one of: 1=1m, 2=3m, 3=10m

    '' radiation sensors are on multiplexer 1
    const MUX1_ENABLED = 0 ' can turn off in the off-season
    const MUX1_QTY = 5      ' 5 radiation sensors to read, plus others...

    '' radiation sensors nameplate data:
    const TLW_SN = ""
    const TLW_MULT = -250
    const ALW_SN = ""
    const ALW_MULT = 215.98
    const ISW_SN = ""
    const ISW_MULT = -101.42
    const RSW_SN = ""
    const RSW_MULT = 112.11
    Const NET_SN = ""
    Const NET_MULT = 14.10  ' output in uV/wm2

    '' evaporation pan pressure transducer:
    const EVAP_MDL = "GlobalWater 4-20mA"
    const EVAP_SN = ""
    const EVAP_MULT = 0.30596 ' -.171
    const EVAP_OFS = -122.385  ' 17.12

    '' thermistor strings reference resistor:
    const Rf_kOhm = 100.3
    const MUX2_QTY = 26

    '' hmp45c AT & RH sensor parameters:
    const RH1_MDL   = "HMP45C"
    const RH1_SN    = "?"
    const RH3_MDL   = "HMP45C"
    const RH3_SN    = "?"
    const RH10_MDL  = "HMP45C"
    const RH10_SN   = "?"

    '' wind sensors:
    const WS1_MDL    = "MetOne 014A"
    const WS1_SN     = "?"
    const WS3_MDL    = "MetOne 014A"
    const WS3_SN     = "?"
    const WS10_MDL   = "RMY 05103"
    const WS10_SN    = "?"

    '' tipping bucket parameters:
    const R_MODL = "TE525MM"
    const R_SN = "?"
    const R_TIP_MM = 0.1  '' mm/tip

    '' radio control parameters
    const RDIO_MDL = "FGR"
    const RDIO_ON_MIN = 0
    const RDIO_OFF_MIN = 15
    const RDIO_V_LIM = 12.2
    const RDIO_HOLD_V_LIM = 12.4

    '' LOCATION PARAMETERS for site (e.g., used for sun angle calcs)
    const LATITUDE = 69.89
    const LONGITUDE = -148.77
    const TZ_HOURS = -9
  
EndConstTable

'' DECLARE FIXED CONSTANTS

'' HMP45C used at 1m, 3m, 10m
const AT_MULT = 0.1 ' degC/mV
Const AT_OFS  = -40 ' degC
Const AT_MULT_10m = 0.14 ' degC/mV
Const AT_OFS_10m  = -80 ' degC
Const ATB_MULT = 1.0
Const ATB_OFS  = 0
const RH_MULT = 0.1 ' %RH/mV
const RH_OFS  = 0   ' %RH

'' MET1 anemometer used at 1m, 3m
const MET1_MULT = 0.0133   ' mps/pulse
const MET1_OFS  = 0.447 ' mps

'' RMYoung wind set used at 10m
const RMY_MULT = 0.098  ' mps/cycle
const RMY_OFS  = 0.01   ' mps (non-zero for WD vector calc if WS=0)
const WD_MULT = 355     ' degrees/fullscale
const WD_OFS  = 0


const MUX_CLOCK_PULSE = 10000 ' 10 msec to step multiplexer relays 

'' DATALOGGER WIRING: channels and ports assignments:

'' SE & DIFF channels:
const AT_1m_SE = 1  ' sensors grouped for use of measurement repetitions
Const AT_3m_SE = 2
Const AT_10m_SE = 3
Const ATB__SE = 8
Const ATB__EX = VX3
const RH_1m_SE = 4
const RH_3m_SE = 5
const RH_10m_SE = 6
const WD_10m_SE = 7

const SLR_PNL_VDIV_SE = 9

Const EVAP_SE = 8
Const NET_DE = 6
const MUX2_SE = 13
const MUX1_DE = 8    ' MUX1 differential input
const MUX1_SEH = 15  ' MUX1 SE inputs are same as DE H & L channels
const MUX1_SEL = 16

'' I/O ports:
const MUX1_ENABLE_PORT = 1   
CONST MUX2_ENABLE_PORT = 2   
CONST MUX_CLOCK_PORT = 3   
CONST HMP_PWR_PORT = 4   
Const SR50_PORT = 5   

'' NOTE: use control+10 for use pulse counting
const WS_3m_PORT = 16    '' cp6
const RAIN_TIP_PORT = 17 '' cp7

const SW12V_USED_FOR = "RF450 radio"

'' pulse counters:
const WS_1m_PULSE = 1
const WS_10m_PULSE = 2

' settings for tag() array, used to log values at startup;
'   public vars tagPtr, biggestTag will show actual use:
const STARTUP_TAGS_COUNT = 46
const STARTUP_TAGS_STR_SIZE = 24

const TYPE_AC = 1
const TYPE_SWITCH = 2
const TYPE_COUNT = 0
const TYPE_FREQ = 1


'' DECLARE PUBLIC VARIABLES (visible via public table)

public diagnostic(3)
alias diagnostic(1) = logger_battery
alias diagnostic(2) = solar_panel
alias diagnostic(3) = panelT
units logger_battery = Volts
units solar_panel = Volts
units PanelT = degC

Public AT(3)
Alias AT(1) = AT_1m
Alias AT(2) = AT_3m
Alias AT(3) = AT_10m
Public ATB_1m
Units ATB_1m = degC

public RH(3)
Units RH() = %
Alias RH(1) = RH_1m
alias RH(2) = RH_3m
alias RH(3) = RH_10m

public DP(3)
Units DP() = degC
Alias DP(1) = DP_1m
alias DP(2) = DP_3m
alias DP(3) = DP_10m

public WS(3)
Units WS() = mps
Alias WS(1) = WS_1m
alias WS(2) = WS_3m
alias WS(3) = WS_10m

public WD
alias WD = WD_10m
units WD = deg

public rain
units rain = mm

public rad_mV(MUX1_QTY)
alias rad_mV(1) = TLW_mV
alias rad_mV(2) = ALW_mV
alias rad_mV(3) = ISW_mV
alias rad_mV(4) = RSW_mV
Alias rad_mV(5) = Net_rad

Public radiation(MUX1_QTY)
Units radiation() = wpm2
Alias radiation(1) = TLW
Alias radiation(2) = ALW
Alias radiation(3) = ISW
Alias radiation(4) = RSW
'alias radiation(5) = Net

Public tc(2)
Units tc() = degC
'Units tc(2) = degC
Alias tc(1) = TLW_dome
Alias tc(2) = ALW_dome

public evap_pan
units evap_pan = cm

Public sr50(4)
Units sr50() = cm
Alias sr50(1) = sr50_raw
alias sr50(2) = sr50_quality
alias sr50(3) = sr50_down
alias sr50(4) = sr50_snowdepth
public sr50_AT ' to be copied from another value!

public thermistor(MUX2_QTY) ' resistance of each sensor, 
units  thermistor = kOhms

'' radio control variables
public RadioOnMin, RadioOffMin
public RadioVlim, RadioHldVlim
Public RadioHold As Boolean
Public RadioForceOn As Boolean

'' startup tags pointer and string size monitor:
public tagPtr, biggestTag


'DECLARE VARIABLES (not visible outside program using 'dim')

dim i
dim ratio
dim tag(STARTUP_TAGS_COUNT) as string * STARTUP_TAGS_STR_SIZE
dim progSig as long

'' note: RealTime(Time()) must be run to update these alias values!
'' Time() array: 1:year, 2:month, 3:day, 4:hour, 5:minute, 6:second,
''             7:microsecond, 8:day_of_week, 9:day_of_year
dim Time(9)
alias Time(4) = thisHour
alias Time(5) = thisMinute
alias Time(9) = thisDayOfYear


'DEFINE DATA TABLES

DataTable(StartupTags,1,10)
    sample(1,progSig,uint2)
    sample(STARTUP_TAGS_COUNT, tag(1), string)
EndTable

DataTable(RainEvnt,Rain,2880)
    sample(1,progSig,uint2)
    sample(1,Rain,fp2)
EndTable

DataTable(Met,1,-1)
    dataInterval(0,60,Min,0)
    Sample(1,progSig,uint2)
    Sample(3,AT,fp2)
    sample(3,RH,fp2)
    sample(3,DP,fp2)
    sample(3,WS,fp2)
    Sample(1,WD,fp2)
    Sample(2,sr50_down,fp2)
    Sample(1,Net_rad,fp2)
    Sample(1,ATB_1m,fp2)
	  'average(1,evap_pan,fp2,0)
EndTable

DataTable(Met_avg,1,-1)
    DataInterval(0,60,Min,0)
    Average(3,AT,fp2,0)
    Average(3,RH,fp2,0)
    Average(3,DP,fp2,0)
    Average(2,WS,fp2,0)
    WindVector(1,WS_10m,WD_10m,fp2,0,0,0,0)
    Totalize(1,rain,fp2,0)
    Average(2,sr50_down,fp2,0)
    Average(1,Net_rad,fp2,0)
    Average(1,ATB_1m,fp2,0)
EndTable

DataTable(Met_min,1,-1)
    dataInterval(0,60,Min,0)
    Minimum(3,AT,fp2,0,0)
    minimum(3,RH,fp2,0,0)
    minimum(3,DP,fp2,0,0)
    minimum(3,WS,fp2,0,0)
    Minimum(2,sr50_down,fp2,0,0)
    Minimum(1,ATB_1m,fp2,0,0)
EndTable

DataTable(Met_max,1,-1)
    dataInterval(0,60,Min,0)
    Maximum(3,AT,fp2,0,0)
    maximum(3,RH,fp2,0,0)
    maximum(3,DP,fp2,0,0)
    maximum(3,WS,fp2,0,0)
    Maximum(2,sr50_down,fp2,0,0)
    maximum(1,ATB_1m,fp2,0,0)
EndTable

DataTable(Logger,1,-1)
    dataInterval(0,60,Min,0)
    sample(1,progSig,uint2)
    sample(3,diagnostic,fp2)
    average(3,diagnostic,fp2,0)
EndTable

DataTable(SR50,1,-1)
    dataInterval(0,5,Min,0)
    sample(1,progSig,uint2)
    sample(1,sr50_AT,fp2)
    sample(4,sr50,fp2)
EndTable

DataTable(radiation,MUX1_ENABLED,-1)
    dataInterval(0,60,Min,0)
    sample(1,progSig,uint2)
    average(MUX1_QTY,radiation,fp2,0)
    average(2,tc,fp2,0)
EndTable

DataTable(thermistors,1,-1)
    dataInterval(0,60,Min,0)
    sample(1,progSig,uint2)
    average(MUX2_QTY,thermistor,fp2,0)
EndTable


' Subroutines

Sub RadioControl( onMin, offMin, batV, Vlimit, holdVlimit )
    '' turn SW12() on & off at specified minutes into interval
    '' battery voltage must be above Vlimit to turn radio on
    '' global RadioHold flag prevents turning radio OFF if set
    '' RadioHold flag is cleared if battery < holdVlimit
    '' special case: radio turns on at noon regardless of V limits
    '' subroutine should normally run every minute

    if RadioForceOn Then SW12(1)

    if thisMinute = onMin Then
        if batV >= Vlimit Then SW12(1)
        if RadioHold AND batV < holdVlimit Then RadioHold = 0
    endIf

    if IfTime(12, 24, hr) Then SW12(1)

    if thisMinute = offMin Then
        if NOT RadioHold Then SW12(0)
    endIf

EndSub ' RadioControl

sub pushTag(name as string * STARTUP_TAGS_STR_SIZE, value as string)
    tagPtr = tagPtr + 1
    if tagPtr > STARTUP_TAGS_COUNT then
        i = tagPtr - STARTUP_TAGS_COUNT
        ' overrun count includes overwritten last entry:
        tag(STARTUP_TAGS_COUNT) = "TAGS_NOT_SAVED=" & i + 1
    else
        tag(tagPtr) = name & "=" & value
        i = len(tag(tagPtr))
        ' see this public variable to optimize tags string size:
        if i > biggestTag then biggestTag = i
    endif
endSub ' pushTag

sub saveStartupTags()
    call pushTag( "SR50_MDL", SR50_MDL )
    call pushTag( "SR50_SN", SR50_SN )
    call pushTag( "SR50_HT_CM", SR50_HT_CM )
    call pushTag( "SR50_AT_SELECT", SR50_AT_SELECT )
    call pushTag( "MUX1_ENABLED", MUX1_ENABLED )
    call pushTag( "MUX1_QTY", MUX1_QTY )
    call pushTag( "TLW_SN", TLW_SN )
    call pushTag( "TLW_MULT", TLW_MULT )
    call pushTag( "ALW_SN", ALW_SN )
    call pushTag( "ALW_MULT", ALW_MULT )
    call pushTag( "ISW_SN", ISW_SN )
    call pushTag( "ISW_MULT", ISW_MULT )
    call pushTag( "RSW_SN", RSW_SN )
    call pushTag( "RSW_MULT", RSW_MULT )
    Call pushTag( "NET_SN", NET_SN )
    Call pushTag( "NET_MULT", NET_MULT )
    call pushTag( "EVAP_MDL", EVAP_MDL )
    call pushTag( "EVAP_SN", EVAP_SN )
    call pushTag( "EVAP_MULT", EVAP_MULT )
    call pushTag( "EVAP_OFS", EVAP_OFS )
    call pushTag( "Rf_kOhm", Rf_kOhm )
    call pushTag( "MUX2_QTY", MUX2_QTY )
    call pushTag( "RH1_MDL", RH1_MDL )
    call pushTag( "RH1_SN", RH1_SN )
    call pushTag( "RH3_MDL", RH3_MDL )
    call pushTag( "RH3_SN", RH3_SN )
    call pushTag( "RH10_MDL", RH10_MDL )
    call pushTag( "RH10_SN", RH10_SN )
    call pushTag( "WS1_MDL", WS1_MDL )
    call pushTag( "WS1_SN", WS1_SN )
    call pushTag( "WS3_MDL", WS3_MDL )
    call pushTag( "WS3_SN", WS3_SN )
    call pushTag( "WS10_MDL", WS10_MDL )
    call pushTag( "WS10_SN", WS10_SN )
    call pushTag( "R_MODL", R_MODL )
    call pushTag( "R_SN", R_SN )
    call pushTag( "R_TIP_MM", R_TIP_MM )
    call pushTag( "RDIO_MDL", RDIO_MDL )
    call pushTag( "RDIO_ON_MIN", RDIO_ON_MIN )
    call pushTag( "RDIO_OFF_MIN", RDIO_OFF_MIN )
    call pushTag( "RDIO_V_LIM", RDIO_V_LIM )
    call pushTag( "RDIO_HOLD_V_LIM", RDIO_HOLD_V_LIM )
    call pushTag( "LATITUDE", LATITUDE )
    call pushTag( "LONGITUDE", LONGITUDE )
    call pushTag( "TZ_HOURS", TZ_HOURS )
endSub ' saveStartupTags

sub read_mux1()

    '' MUX1 has 5 radiometers, 2 thermocouples, 1 pressure transducer,
    '' and a 10:1 voltage divider
    
    PortSet(MUX1_ENABLE_PORT, 1)
    
    '' terrestrial long wave radiation
    PulsePort(MUX_CLOCK_PORT ,MUX_CLOCK_PULSE)
    VoltDiff(TLW,1,mv25,MUX1_DE,1,0,250,TLW_MULT,0)
    
    '' atmospheric long wave radiation
    PulsePort(MUX_CLOCK_PORT ,MUX_CLOCK_PULSE)
    VoltDiff(ALW,1,mv25,MUX1_DE,1,0,250,ALW_MULT,0)
    
    '' incident short wave radiation
    PulsePort(MUX_CLOCK_PORT ,MUX_CLOCK_PULSE)
    VoltDiff(ISW,1,mv25,MUX1_DE,1,0,250,ISW_MULT,0)
    
    '' reflected short wave radiation
    PulsePort(MUX_CLOCK_PORT ,MUX_CLOCK_PULSE)
    VoltDiff(RSW,1,mv25,MUX1_DE,1,0,250,RSW_MULT,0)
    
    '' net radiation
    PulsePort(MUX_CLOCK_PORT ,MUX_CLOCK_PULSE)
		VoltDiff(Net_rad,1,mv25,NET_DE,True,0,250,1000/NET_MULT,0)
		  If WS_3m>=5 Then
		   Net_rad=Net_rad*(1+0.021286*(WS_3m-5))
      Else
	     Net_rad=Net_rad
	    EndIf
    
    '' IMPORTANT: thermocouple reference temperature is currently
    '' set to use logger panel temp, but this assumes that the 
    '' multiplexer is in the same enclosure with the same temperature
    '' environment and gradients.
    
    '' measure radiometer thermocouples:
    PulsePort(MUX_CLOCK_PORT ,MUX_CLOCK_PULSE)
    TCDiff (tc(1),1,mV2_5C,MUX1_DE,TypeT,panelT,1,0,250,1.0,0)
    
    PulsePort(MUX_CLOCK_PORT ,MUX_CLOCK_PULSE)
    TCDiff (tc(2),1,mV2_5C,MUX1_DE,TypeT,panelT,1,0,250,1.0,0)
    
    PortSet(MUX1_ENABLE_PORT, 0)

EndSub ' read_mux1


BeginProg

    '' NOTE: STARTUP CODE RUNS ONCE BEFORE SCANS BEGIN!
    progSig = status.ProgSignature
    Call saveStartupTags()
    CallTable StartupTags

    '' turn radio and RadioHold on when the program starts:
    SW12(1)
    RadioHold = 1

    '' copy initial settings from contstants into variables:
    RadioOnMin = RDIO_ON_MIN
    RadioOffMin = RDIO_OFF_MIN
    RadioVlim = RDIO_V_LIM
    RadioHldVlim = RDIO_HOLD_V_LIM

    '' END OF STARTUP CODE

    Scan(60,Sec,0,0)     'SCAN PROGRAM AT 60-Second INTERVALS

        RealTime(Time())       '' sets thisMinute and other alias variables

        '' read logger diagnostic values:

        Battery(logger_battery)
        VoltSe(solar_panel,1,mV5000,SLR_PNL_VDIV_SE,1,0,250,0.01,0)
        PanelTemp(panelT,250)

        '' read HMP45C AT & RH sensors:
        ''    turn on relay to power to HMP45C
        ''    delay for sensors startup
        ''    measure AT, RH
        ''    calculate dewpoint
        ''    turn off relay

        PortSet(HMP_PWR_PORT,1 )
        Delay(0,4000,mSec) ' delay for HMP sensor startup
        VoltSe(AT_1m,1,mV2500,AT_1m_SE,1,0,250,AT_MULT,AT_OFS)
        VoltSe(AT_3m,1,mV2500,AT_3m_SE,1,0,250,AT_MULT,AT_OFS)
        VoltSe(AT_10m,1,mV2500,AT_10m_SE,1,0,250,AT_MULT_10m,AT_OFS_10m)
        VoltSe(RH_1m,1,mV2500,RH_1m_SE,1,0,250,RH_MULT,RH_OFS)
        VoltSe(RH_3m,1,mV2500,RH_3m_SE,1,0,250,RH_MULT,RH_OFS)
        VoltSe(RH_10m,1,mV2500,RH_10m_SE,1,0,250,RH_MULT,RH_OFS)
        VoltSE(evap_pan,1,mv2500,EVAP_SE,1,0,250,EVAP_MULT,EVAP_OFS)
        dewPoint(DP(1),AT(1),RH(1))
        dewPoint(DP(2),AT(2),RH(2))
        dewPoint(DP(3),AT(3),RH(3))
        PortSet(HMP_PWR_PORT,0)
        
		    '107 Temperature Probe measurement 'atb'
		    Therm107(ATB_1m,1,ATB__SE,ATB__EX,0,250,ATB_MULT,ATB_OFS)

          '' count rain bucket-tip events since last scan:
        PulseCount(rain,1,RAIN_TIP_PORT,TYPE_SWITCH,0,R_TIP_MM,0)

        '' measure wind speed and direction:
        pulseCount(WS(1),1,WS_1m_PULSE, TYPE_SWITCH,TYPE_COUNT,MET1_MULT, MET1_OFS)
        pulseCount(WS(2),1,WS_3m_PORT,  TYPE_SWITCH,TYPE_COUNT,MET1_MULT, MET1_OFS)
        pulseCount(WS(3),1,WS_10m_PULSE,TYPE_AC,TYPE_FREQ,RMY_MULT,RMY_OFS)
        brHalf(WD,1,mV2500,WD_10m_SE,Vx1,1,2500,1,0,250,WD_MULT,WD_OFS)

        '' read SR50 snow depth sensor:
        ''    read distance (m>0) and quality number into sr50_raw(2)
        ''    apply speed-of-sound correction and convert to cm down
        ''    calculate snow depth as sum of distance down and sensor height

        SDI12Recorder(sr50, SR50_PORT, SR50_SDI_ADDR, "M1!", 1, 0)
        sr50_AT = AT(SR50_AT_SELECT)
        sr50_down = sr50 * -100 * SQR((sr50_AT+273.15) / 273.15)
        sr50_snowdepth = sr50_down + SR50_HT_CM

        'NR-LITE2 Net Radiometer (dynamic wind speed correction) measurement 'nr' and 'corr_nr'
		    VoltDiff(Net_rad,1,mv25,NET_DE,True,0,250,1000/NET_MULT,0)
		    If WS_3m>=5 Then
			    Net_rad=Net_rad*(1+0.021286*(WS_3m-5))
		    Else
			     Net_rad=Net_rad
		    EndIf

        callTable RainEvnt
        callTable Logger
        callTable Met
        callTable Met_avg
        callTable Met_min
        callTable Met_max
        callTable SR50

        if MUX1_ENABLED then call read_mux1()

        callTable Radiation

        call RadioControl(RadioOnMin, RadioOffMin, logger_battery, RadioVlim, RadioHldVlim )
    
    nextScan

    slowsequence '' read thermistor strings every 10 minutes:
        scan(600,Sec,0,0)

            '' NOTE: I think we can rely on MUX1 already been finished,
            '' otherwise use of common clock signal might be a problem...

            '' process thermistors:
            ''    turn on multiplexer
            ''    iterate MUX2_QTY times
            ''    read 1 thermistor at each step
            ''    calculate thermistor resistance from ratio and ref R
            ''    turn off multiplexer when done

            PortSet(MUX2_ENABLE_PORT, 1)
            i = 1
            subScan(0,usec,MUX2_QTY)
                PulsePort(MUX_CLOCK_PORT, MUX_CLOCK_PULSE)
                BrHalf (ratio,1,mV2500,MUX2_SE,Vx2,1,2500,1,0,250,1.0,0)
                thermistor(i) = Rf_kOhm * ratio / (1-ratio) 
                i = i + 1
            nextSubScan
            PortSet(MUX2_ENABLE_PORT, 0)
    
            calltable thermistors

        nextscan
    endsequence

EndProg