UTESpac/UTES_struct_setup.m at master · UtahEFD/UTESpac · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
%Prepare to calculate 2nd order structure parameters of temperature, humidity, and
%the temperature humidity cross structure parameter

%Can add to calculate further structure parameters here at a later date.

function [output, raw] = UTES_struct_setup(data, info, output, raw, tableNames, sensorInfo)

if isempty(raw)
    warning('Could not calculate Structure parameters. raw data does not exist')
    return
end

fprintf('\nCalculating Structure Parameters\nC_T^2 C_Tq C_q^2\n');

Rd = 287.058;  % [J/K/kg] Gas constant for air


%%%%%%%%%%%%%
%When Finewire does not exist
%Eventually need a better system than this
fields = fieldnames(sensorInfo);
no_fw = ~any(cellfun(@(x) strcmp(x, 'fw'), fields));
no_irga = ~any(cellfun(@(x) strcmp(x, 'irgaH2O'), fields));
%%%%%%%%%%%%%


%Number of Sonics
NumSonic = size(sensorInfo.Tson, 1);
%Number IRGAs
if no_irga
    NumIRGA = 0;
else
    NumIRGA = size(sensorInfo.irgaH2O, 1);
end

MaxInst = max([NumSonic, NumIRGA]);
MinInst = min([NumSonic, NumIRGA]);


%Iterate through number of instruments (heights & Tables)
cntrIRGA = 1;
for qq=1:MaxInst
    %TableName columns in UTESpac and siteInfo do not match....
    %These lines look for the tablename from siteInfo
    TableNum = sensorInfo.Tson(qq, 1);
    tmp = regexp(info.tableNames, tableNames{TableNum});
    [~, siteInfoTableNum] = max(~cellfun(@isempty, tmp));

    %Calculate number of points per averaging period
    pnts = info.tableScanFrequency(:, siteInfoTableNum)*60*info.avgPer;

    %Create average time vector
    time = raw.t(pnts:pnts:end);

    %Iterate through each aeraging period and calculate strucutre
    %parameters
    cntr1=1;
    cntr2 = pnts;
    for ii=1:length(raw.t)/pnts

        % find reference pressure
        zRef = min(sensorInfo.u(:,3));  % zRef is lowest sonic level
        altitude = info.siteElevation;
        PRef = 101325*(1-2.25577*10^-5*(altitude+zRef))^5.25588; % find pRef from elevation: http://www.engineeringtoolbox.com/air-altitude-pressure-d_462.html


        SonHeight = sensorInfo.Tson(qq, 3);
        %Sonic Variables for each height
        sigma_u = nanstd(raw.uPF(cntr1:cntr2, qq));
        sigma_v = nanstd(raw.vPF(cntr1:cntr2, qq));
        sigma_w = nanstd(raw.wPF(cntr1:cntr2, qq));
        U = nanmean(sqrt(raw.uPF(cntr1:cntr2, qq).^2+raw.vPF(cntr1:cntr2, qq).^2));
        if isfield(sensorInfo, 'P')
            if and(size(sensorInfo.P, 1)==1, sensorInfo.P(1)==siteInfoTableNum)
                P = nanmean(data{TableNum}(cntr1:cntr2, sensorInfo.P(1, 2)).*1000); %average Pressure [Pa]
            else
                P = nanmean(data{TableNum}(cntr1:cntr2, sensorInfo.P(siteInfoTableNum, 2)).*1000); %average Pressure [Pa]
            end

            if isnan(P) || abs((P - PRef)/PRef) > 0.05
                P = PRef;
            end
        else
            P = PRef;
        end

        if ~no_fw
            T = raw.fwT(cntr1:cntr2, qq)+273.15; %Temperature [K]

            %Calculate temperature structure parameter
            [~, Ct2(ii)] = ...
                crossStruct(T, T, 20, 'spatial', sigma_u, sigma_v, sigma_w, U);
            if Ct2(ii)<0
                Ct2(ii) = nan;
            end
        else
            T = nan;
        end
        TSon = raw.sonTs(cntr1:cntr2, qq)+273.15; %Sonic Temperature [k]


        %Calculate temperature structure parameter with Sonic temperature
        [~, Ctson2(ii)] = ...
            crossStruct(TSon, TSon, 20, 'spatial', sigma_u, sigma_v, sigma_w, U);
        if Ctson2(ii)<0
            Ctson2(ii) = nan;
        end


        if cntrIRGA<=MinInst
            %look for matching IRGA for Sonic height
            IRGAHeight = find(sensorInfo.irgaH2O(:, 3)==SonHeight);

            if ~isempty(IRGAHeight)
                %Abslute humidity [kg m^-3]
                rhov = raw.rhov(cntr1:cntr2, IRGAHeight)./1000;

                %Calculate mean temperature for period
                T_mean = nanmean(T);
                if isnan(T_mean)
                    T_mean = nanmean(TSon);
                end

                %Specific humidity [kg/kg]
                q = (rhov./(P./(Rd*T_mean)));

                %Calculate temperature humidity cross structure parameter


                %Calculate humidity structure parameter
                [~, Cq2(ii)] = ...
                    crossStruct(q, q, 20, 'spatial', sigma_u, sigma_v, sigma_w, U);
                if Cq2(ii)<0
                    Cq2(ii) = nan;
                end

                %Calculate temperature humidity correlation
                if ~no_fw
                    [~, Ctq(ii)] = ...
                        crossStruct(T, q, 20, 'spatial', sigma_u, sigma_v, sigma_w, U);

                    r_tq(ii) = Ctq(ii)./sqrt(Ct2(ii).*Cq2(ii));
                end

                %Calculate temperature humidity correlation TSonic
                [~, Ctsonq(ii)] = ...
                    crossStruct(TSon, q, 20, 'spatial', sigma_u, sigma_v, sigma_w, U);

                %Calculate temperature humidity correlation
                r_tsonq(ii) = Ctsonq(ii)./sqrt(Ctson2(ii).*Cq2(ii));

                IRGAFlag = 1;
            else
                IRGAFlag = 0;
            end
        else
            IRGAFlag = 0;
        end
        %update counters
        cntr1 = cntr1 + pnts;
        cntr2 = cntr2 + pnts;
    end

    %Save to output structure
    if qq==1
        if IRGAFlag
            if no_fw
                %If IRGA Exists @ no finewire
                output.StructParam = [time, Ctson2', Ctsonq', Cq2', r_tsonq'];
                output.StructParamHeader = {'Timestamp', ...
                    [num2str(sensorInfo.Tson(qq, end)), 'm CtSon2: K^2 m^-2/3'],...
                    [num2str(sensorInfo.irgaH2O(cntrIRGA, end)), 'm CtSonq: K kg kg^-1 m^-2/3'],...
                    [num2str(sensorInfo.irgaH2O(cntrIRGA, end)), 'm Cq2: kg^2 kg^-2 m^-2/3'],...
                    [num2str(sensorInfo.irgaH2O(cntrIRGA, end)), 'm r_tSonq'];...
                    [],num2str(sensorInfo.Tson(qq, end)),...
                    num2str(sensorInfo.Tson(qq, end)),...
                    num2str(sensorInfo.irgaH2O(cntrIRGA, end)),...
                    num2str(sensorInfo.irgaH2O(cntrIRGA, end))};
                cntrIRGA = cntrIRGA+1;
            else
                %If IRGA Exists and FW exists
                output.StructParam = [time, Ct2', Ctson2', Ctq', Ctsonq', Cq2', r_tq', r_tsonq'];
                output.StructParamHeader = {'Timestamp', ...
                    [num2str(sensorInfo.Tson(qq, end)), 'm Ctfw2: K^2 m^-2/3'],...
                    [num2str(sensorInfo.Tson(qq, end)), 'm CtSon2: K^2 m^-2/3'],...
                    [num2str(sensorInfo.irgaH2O(cntrIRGA, end)), 'm Ctfwq: K kg kg^-1 m^-2/3'],...
                    [num2str(sensorInfo.irgaH2O(cntrIRGA, end)), 'm CtSonq: K kg kg^-1 m^-2/3'],...
                    [num2str(sensorInfo.irgaH2O(cntrIRGA, end)), 'm Cq2: kg^2 kg^-2 m^-2/3'],...
                    [num2str(sensorInfo.irgaH2O(cntrIRGA, end)), 'm r_tfwq'],...
                    [num2str(sensorInfo.irgaH2O(cntrIRGA, end)), 'm r_tSonq'];...
                    [],num2str(sensorInfo.Tson(qq, end)),...
                    num2str(sensorInfo.Tson(qq, end)),...
                    num2str(sensorInfo.irgaH2O(cntrIRGA, end)),...
                    num2str(sensorInfo.irgaH2O(cntrIRGA, end)),...
                    num2str(sensorInfo.irgaH2O(cntrIRGA, end)),...
                    num2str(sensorInfo.irgaH2O(cntrIRGA, end)),...
                    num2str(sensorInfo.irgaH2O(cntrIRGA, end))};
                cntrIRGA = cntrIRGA+1;
            end
        else
            %If only Sonics
            output.StructParam = [time, Ct2', Ctson2'];
            output.StructParamHeader = {'Timestamp', ...
                [num2str(sensorInfo.Tson(qq, end)), 'm Ctfw2: K^2 m^-2/3'],...
                [num2str(sensorInfo.Tson(qq, end)), 'm CtSon2: K^2 m^-2/3'];...
                [],num2str(sensorInfo.Tson(qq, end)),...
                num2str(sensorInfo.Tson(qq, end))};
        end
    else
        if IRGAFlag
            if no_fw
                %IRGA Exists
                output.StructParam = [output.StructParam, Ctson2', Ctsonq', Cq2', r_tsonq'];
                output.StructParamHeader = [output.StructParamHeader,...
                    {[num2str(sensorInfo.Tson(qq, end)), 'm CtSon2: K^2 m^-2/3'],...
                    [num2str(sensorInfo.irgaH2O(cntrIRGA, end)), 'm CtSonq: K kg kg^-1 m^-2/3'],...
                    [num2str(sensorInfo.irgaH2O(cntrIRGA, end)), 'm Cq2: kg^2 kg^-2 m^-2/3'],...
                    [num2str(sensorInfo.irgaH2O(cntrIRGA, end)), 'm r_tSonq'];...
                    num2str(sensorInfo.Tson(qq, end)),...
                    num2str(sensorInfo.Tson(qq, end)),...
                    num2str(sensorInfo.irgaH2O(cntrIRGA, end)),...
                    num2str(sensorInfo.irgaH2O(cntrIRGA, end))}];
                cntrIRGA = cntrIRGA+1;
            else
                %IRGA Exists
                output.StructParam = [output.StructParam, Ct2', Ctson2', Ctq', Ctsonq', Cq2', r_tq',r_tsonq'];
                output.StructParamHeader = [output.StructParamHeader,...
                    {[num2str(sensorInfo.Tson(qq, end)), 'm Ctfw2: K^2 m^-2/3'],...
                    [num2str(sensorInfo.Tson(qq, end)), 'm CtSon2: K^2 m^-2/3'],...
                    [num2str(sensorInfo.irgaH2O(cntrIRGA, end)), 'm Ctfwq: K kg kg^-1 m^-2/3'],...
                    [num2str(sensorInfo.irgaH2O(cntrIRGA, end)), 'm CtSonq: K kg kg^-1 m^-2/3'],...
                    [num2str(sensorInfo.irgaH2O(cntrIRGA, end)), 'm Cq2: kg^2 kg^-2 m^-2/3'],...
                    [num2str(sensorInfo.irgaH2O(cntrIRGA, end)), 'm r_tfwq'],...
                    [num2str(sensorInfo.irgaH2O(cntrIRGA, end)), 'm r_tSonq'];...
                    num2str(sensorInfo.Tson(qq, end)),...
                    num2str(sensorInfo.Tson(qq, end)),...
                    num2str(sensorInfo.irgaH2O(cntrIRGA, end)),...
                    num2str(sensorInfo.irgaH2O(cntrIRGA, end)),...
                    num2str(sensorInfo.irgaH2O(cntrIRGA, end)),...
                    num2str(sensorInfo.irgaH2O(cntrIRGA, end)),...
                    num2str(sensorInfo.irgaH2O(cntrIRGA, end))}];
                cntrIRGA = cntrIRGA+1;
            end
        else
            %Only Sonics
            output.StructParam = [output.StructParam, Ct2', Ctson2'];
            output.StructParamHeader = [output.StructParamHeader, ...
                {[num2str(sensorInfo.Tson(qq, end)), 'm Ctfw2: K^2 m^-2/3'],...
                [num2str(sensorInfo.Tson(qq, end)), 'm CtSon2: K^2 m^-2/3'];...
                num2str(sensorInfo.Tson(qq, end)),...
                num2str(sensorInfo.Tson(qq, end))}];
        end
    end

    clearvars C* r_*

end