-
Notifications
You must be signed in to change notification settings - Fork 0
/
larpkg.ipf
6354 lines (5704 loc) · 228 KB
/
larpkg.ipf
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
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
//=====================================================================================
// LAR Library
//
// Patrick O'Keeffe Laboratory for Atmospheric Research
// pokeeffe@wsu.edu Dept. of Civil & Environmental Engineering
// 509.335.7246 Washington State University
//
// Licensed under the MIT License, available at http://opensource.org/licenses/MIT
//
//=====================================================================================
#pragma rtGlobals = 1 // Use modern global access method.
#pragma IgorVersion = 6.20
#pragma version = 0.20130313 // this field last changed on, YYYYMMDD
#If ( 1 ) // 1 to compile; 0 to disable procedure
//=====================================================================================
// INCLUDES
//=====================================================================================
#include <Remove Points>
//#include <New Polar Graphs>
//#include <Rose Plot>
//#include <WaveSelectorWidget> ??
//#include <Time-Frequency> ??
//====================================================================================
// CONSTANTS
//=====================================================================================
StrConstant ksFileFilter = "Plain-text data (.txt .dat .csv .log):.txt,.dat,.csv,.log;Archived data (.zip .7z .rar):.zip,.7z,.rar;All files (*.*):.*;"
// molecular weights
Constant kMW_dryair = 28.9645 // g / mol dry air
//Constant kMW_ch4 =
Constant kMW_co2 = 44.01 // g / mol carbon dioxide
//Constant kMW_co =
Constant kMW_h2o = 18.01528 // g / mol water
//Constant kMW_n2o =
// gas constants
Constant kRu = 8.3144621e-2 // (mb m^3) / (mol K) universal gas constant
Constant kR_dryair = 2.87057e-3 // (mb m^3) / (g K) specific gas constant for dry air
Constant kR_h2o = 4.61495e-3 // (mb m^3) / (g K) specific gas constant for water vapor
// specific heats
Constant kCp_dryair = 1.0057e-2 // (mb m^3) / (g K) specific heat capacity at constant pressure for dry air
//Constant kCpv = 1875 // ± 25 J / (kg K) specific heat capacity at constant pressure for water vapor ????
// other constants
Constant kVonKarman = 0.4 // dimensionless
Constant kGravity = 9.8 // m / s^2
// as defined for use in structures
Constant MAX_OBJ_NAME = 31 // names
Constant MAX_WIN_PATH = 199 // windows
Constant MAX_UNITS = 49 // axis units
Constant MAXCMDLEN = 400 // string values
//=====================================================================================
// PUBLIC FUNCTIONS
//=====================================================================================
// returns <wrefs> with new references <addref> at point <point>
Function/WAVE AddWaveRef( addref, wrefs, beforePoint )
wave addref // reference to add to <wrefs>
wave/WAVE wrefs // wave of references to add to
variable beforePoint // inserted at this point, existing values moved backward
// 0 = front numpnts(wrefs) = back
InsertPoints beforePoint, 1, wrefs
wrefs[beforePoint] = addref
return wrefs
end
// Returns wave of references to all data folders in current data folder (except Packages, if in root:)
//
// 2011.09 initial release
Function/WAVE AllDataFoldersHere(sortBy)
variable sortBy
string found = DataFolderDir(1)
found = RemoveFromList("Packages", ReplaceString(",", found[8, strlen(found)-2], ";") )
if ( sortBy+1 ) // true except for -1
found = SortList( found, ";", sortBy)
endif
Make/FREE/DF/N=(ItemsInList(found)) reflist
variable i
for (i=0; i<ItemsInList(found); i+=1)
reflist[i] = $StringFromList(i, found)
endfor
return reflist // not threadsafe
End
// Returns wave containing references to all waves in current data folder
//
// 2011.09 initial release
Function/WAVE AllWavesHere(sortBy)
variable sortBy
string found = WaveList("*", ";", "")
if ( sortBy+1 ) // true except for -1
found = SortList( found, ";", sortBy )
endif
Make/FREE/WAVE/N=(ItemsInList(found)) refw
variable i
for (i=0; i<ItemsInList(found); i+=1)
refw[i] = $StringFromList(i, found)
endfor
return refw // not threadsafe
End
// returns ambient temp, Celcius, derived from sonic temp
//
// 2011.11.10 now works in Celcius
// 2011.11.09 written
ThreadSafe Function AmbientTemp( Ts, Q_ )
variable Ts // abs. sonic temp. Celcius
variable Q_ // specific humidity dimensionless
return ( (Ts+273.15) / (1+ 0.51*Q_) ) - 273.15
End
// writes all the waves in <wrefs> to a .csv for each interval of size <interval> in <tstamp>
//
// 2011.09 initial release
Function BaleWaves( tstamp, wrefs, interval, aligned, formatTableName, destNameMask, destPath, overwrite)
wave/D tstamp // time to bale according to
wave/WAVE wrefs // waves to affect, in desired order of output L->R
variable interval // bale interval in seconds
variable aligned // nonzero to start at some multiple of the interval
string formatTableName // name of existing table with desired formatting
string destNameMask // output file name format; use field codes from StringFromMaskedVar
string destPath // full path to output directory or "" for prompt
variable overwrite // nonzero to overwrite existing files, default: get prompted to save as...
variable oi, wn, mgc, ints // output index, wave num, modifygraph count, intervals
string tmpTN, mglist, outname // temp table name, modifygraph list, output file name
DoWindow $formatTableName // not threadsafe
If ( !V_flag )
print ("BaleWaves: could not find table <"+formatTableName+"> - aborting")
return -1
endif
If (strlen(destPath))
NewPath/Q/O/M=("Select output directory") balewoutdir, destPath
else
NewPath/Q/O/M=("Select output directory") balewoutdir
endif
wave/D bounds = IntervalBoundaries( tstamp, interval, aligned)
ints = DimSize(bounds,0)
DFREF savDFR = GetDataFolderDFR()
SetDataFolder NewDataFolderX("root:Packages:art:BaleWavesTemp")
mglist = ReplaceString(num2char(13),WinRecreation(formatTableName, 1),";")
for (oi=0; oi<ints; oi+=1)
Edit/HIDE=1
tmpTN = S_name // returned by Edit
for (wn=0; wn<numpnts(wrefs); wn+=1)
wave now = wrefs[wn]
Duplicate/O/R=[ bounds[oi][%lo], bounds[oi][%hi] ] now, $NameOfWave(now)
AppendToTable $NameOfWave(now)
endfor
// mglist = used to be here... haven't tested since moved outside loop!
for (mgc=3; mgc<(ItemsInList(mglist)-2); mgc+=1)
Execute StringFromList(mgc, mglist)
endfor
outname = StringFromMaskedVar( destNameMask, tstamp[ bounds[oi][%lo] ] )
If (overwrite)
SaveTableCopy/O/W=$tmpTN/T=2/P=balewoutdir as outname
else
SaveTableCopy/W=$tmpTN/T=2/P=balewoutdir as outname
endif
If ( V_flag )
DoWindow/K $tmpTN
print "BaleWaves: Detected error condition or user cancel - aborting"
break
endif
DoWindow/K $tmpTN
endfor
KillDataFolder :
SetDataFolder savDFR
printf "BaleWaves: finished writing %d files ", ints
End
// Rounds value <inVal> to specified digits column <place> using bankers' rounding; ie: if the remainder
// is one-half (0.5), rounding proceeds towards the nearest even integer rather than upwards to avoid an
// upward statistical bias over large data sets.
//
// 2011.10.27 *bug fix* added second abs() around statement compared to maxdev
// 2010.old initial release
Function BankerRound( inVal, place, [toOdd] )
variable inVal // floating point
variable place // digits column to round to (0=1's place, 3= 10^3 =1000s, -3= 10^-3 =1/1000ths)
variable toOdd // non-zero for round-to-odd
variable maxdev = ( 0.00001 < 10^(place-5) ? 0.00001 : 10^(place-5) )
// R ~= 1/2? chose 5+ places smaller than precision, at least within 1/100,000th of 1/2
variable R = mod(inVal, 10^place)/10^place // calc normalized remainder R (range 0 to +/-1)
if ( abs(abs(R) - 0.5) <= maxdev ) // if very close to +/-0.5
// make inVal (+) integer, modulo by 2 to determine inVal as even/odd (0=even, 1=odd)
// evaluate toOdd boolean expression (0=round-to-even, 1=round-to-odd)
// if inVal=even & round-to-even OR inVal=odd & round-to-odd, nudge R closer to zero
// if inVal=odd & round-to-even OR inVal=even & round-to-odd, nudge R away from zero
R += ( mod(trunc(abs(inVal)),2)==(toOdd!=0) ) ? -0.2*sign(R) : 0.2*sign(R)
endif
inVal /= 10^place // move decimal to target place
inVal = ( abs(R) < 0.5 ) ? trunc(inVal) : round(inVal) // if |R|<1/2, round to 0, otherwise round away 0
inVal *= 10^place // move decimal back to starting place & return
return inVal // return rounded inVal
End
// Returns a string with the bits in <var> written MSB->LSB as 4bit blocks
// arbitrarily limited to 32 bits. change variable <upperBound> if desired
//
// 2013.05.23 major API behavior & changes
// - arg `howMany` has been removed; function automatically determines length of
// resulting string now
// - new arg `minLen` can be used to get leading zeros if desired
// - new arg `wordsize` determines grouping within string (1110 0100 vs 11 10 01 00)
// defaults to four
// - arg `maxlen` operates generally as expected; also, it overrides minlen (!) if max<min
// 2011.09 initial release
Function/S BitString( num [, wordsize, minLen, maxLen] )
variable num, wordsize, minLen, maxLen
variable i
string out = ""
num = trunc(num)
wordsize = ( ParamIsDefault(wordsize) ? 4 : wordsize )
minLen = Max((ParamIsDefault(minLen) ? 0 : minLen), ceil(log(num)/log(2)))
maxLen = ( ParamIsDefault(maxLen) ? minLen : maxLen )
for ( i=Limit(minLen, 1, maxLen)-1; i>=0; i-=1 )
sprintf out, "%s%d", out, ( (num & (2^i)) != 0 )
if ( !Mod( i, wordsize ) )
out += " "
endif
endfor
return out // not threadsafe
End
// Returns numeric interpretation of identifiable cardinal directions (NW, S, SSE) or NAN
//
// 2011.09 initial release
Function Cardinal2D( inStr )
string inStr // not threadsafe
variable outVal
inStr = UpperStr(UnPadString(inStr, 0x20)) // remove trailing spaces, convert uppercase
strswitch(inStr)
case "N":
return 0
case "NNE":
return 22.5
case "NE":
return 45
case "ENE":
return 67.5
case "E":
return 90
case "ESE":
return 112.5
case "SE":
return 135
case "SSE":
return 157.5
case "S":
return 180
case "SSW":
return 202.5
case "SW":
return 225
case "WSW":
return 247.5
case "W":
return 270
case "WNW":
return 292.5
case "NW":
return 315
case "NNW":
return 337.5
default:
return NAN
endswitch
end
// Returns <var> with bit # <bit> set to 0
//
// derived directly from the demo under Using Bitwise Operators
//
// 2011.09 initial release
Function ClearBit( var, bit )
variable var, bit
return ( trunc(var) & ~(2^bit) )
End
// TODO
// consider re-implementing the X-scale check
// add a wtype filter to pass-thru to WaveList
//
// concatenates waves of the same name from datafolders in <dfrlist>, placing them in <destDFR>
//
// 2012.02.28 added feedback to differentiate between intial and subsequent additions to concat list;
// added alert for empty (skipped) folders
// 2011.11.02 fixed bug where folders were not killed with <kill=1> and added an empty-check
// 2011.10.28 changed output destination from DFREF to string path
// 2010.old derived from old function
Function ConcatAcrossDFRs( DFRlist, destPath, overwrite, kill, [wfilter] )
wave/DF DFRlist // ordered list of data folders to concatenate from
string destPath // string path to output folder, passed internally to NewDataFolderX
// ""
variable overwrite // nonzero to overwrite existing data in output folder
variable kill // 0 do not kill sources
// 1 kill sources after concatenation (memory conserving)
// 2 kill sources after finishing all concatenations (for paranoid androids)
string wfilter // optional filter passed to wavelist(); defaults to "*"
wfilter = SelectString(ParamIsDefault(wfilter), wfilter, "*")
DFREF sav0 = GetDataFolderDFR()
DFREF dest = NewDataFolderX( destPath )
variable i, j, nwaves, nlists, windex, endpnt, initialAdditions
string wfound, wname, thislist
Make/T/FREE/N=0 concatlist
for (i=0; i<numpnts(DFRlist); i+=1)
SetDataFolder DFRlist[i]
wfound = WaveList(wfilter,";","")
nwaves = ItemsInList(wfound)
if ( !nwaves )
print "* Warning: folder <"+GetDataFolder(1, DFRlist[i])+"> contained no waves - skipping"
continue // skip folder if no waves inside
else
print "Adding folder <"+GetDataFolder(1, DFRlist[i])+"> to source list"
endif
if ( numpnts(concatlist) )
initialAdditions = 1
endif
for (j=0; j<nwaves; j+=1)
wname = StringFromList(j, wfound)
windex = FindDimLabel(concatlist, 0, wname)
If (windex < -1)
endpnt = numpnts(concatlist)
InsertPoints endpnt, 1, concatlist
SetDimLabel 0, endpnt, $wname, concatlist
concatlist[endpnt] = GetDataFolder(1)+wname+";"
If ( initialAdditions )
print "\t* Warning: added previously unfound wave to list: "+wname
else
print "\tAdded wave to concatenation list: "+wname
endif
else
concatlist[windex] += GetDataFolder(1)+wname+";"
endif
endfor
endfor
SetDataFolder sav0
for (i=0; i<numpnts(concatlist); i+=1)
thislist = concatlist[i]
// // EXCERPTED FROM OLD CODE - RELIABILITY SUSPECT
// if ( checkTime ) // if supposed to verify continuous timing
// timeOK = 1 // assume true
// for (j=1; j<ItemsInList(thisList); j+=1) // for each pair of waves in the list
// wave A = $StringFromList(j-1, thisList) // build ref to earlier
// wave B = $StringFromList(j, thisList) // build ref to later
// if ( (rightx(A)!=leftx(B)) || (deltax(A)!=deltax(B)) ) // if discontinuous or diff timestep
// // lar_print(2, "*\tWaves "+GetWavesDataFolder(A,2)+" and "+GetWavesDataFolder(B,2)+" had discontinuous times and were not concatenated")
// timeOK = 0 // display msg, set false
// break // leave j-for loop
// endif
// endfor
// if ( !timeOK ) // if a time discontinuity was found
// // lar_print(0, ">\tWaves of the name "+NameOfWave(A)+" were not concatenated because of time discontinuities.")
// continue // skip remainder of i-for loop
// endif
// endif
// // END EXCERPT
wname = GetDimLabel(concatlist, 0, i)
If (overwrite && kill == 1)
Concatenate/NP/KILL/O thislist, dest:$wname
elseif (overwrite && kill != 1)
Concatenate/NP/O thislist, dest:$wname
elseif (!overwrite && kill ==1)
Concatenate/NP/KILL thislist, dest:$wname
else
Concatenate/NP thislist, dest:$wname
endif
endfor
If (kill)
for (i=0; i<numpnts(DFRlist); i+=1)
DFREF dfr = DFRlist[i]
If (CountObjectsDFR(dfr,1) || CountObjectsDFR(dfr,2) || CountObjectsDFR(dfr,3) || CountObjectsDFR(dfr,4) )
print "ConcatAcrossDFRs: could not kill data folder <"+GetDataFolder(0,dfr)+"> because it still contains objects"
else
KillDataFolder/Z DFRlist[i]
endif
endfor
endif
return 0
End
// Returns total number of NANs in wave
//
// 2011.10.11 added Limit()s to point boundaries
// 2011.09 initial release
Function CountNans( wname, [p1, p2] )
wave wname // wave to examine
variable p1, p2 // optional point boundaries, inclusive
variable i, tot
p1 = Limit(p1, 0, p1)
p2 = Limit( (ParamIsDefault(p2) ? numpnts(wname)-1 : p2), p1, numpnts(wname)-1 )
for (i=p1; i<=p2; i+=1)
tot += (numtype(wname[i])==2)
endfor
return tot
End
// Returns covariance of wx and wy as Cov(wx,wy) = mean(wx*wy) - mean(wx)*mean(wy)
// or NAN if any points were NAN or waves not same length
//
// 2011.11.21 changed numpnts() to DimSize(); added multidimensional check
// 2011.10.11 added limits to point boundaries
// 2011.10 added point boundaries
// 2011.09 initial release
Function Cov( wx, wy, [p1, p2])
wave wx, wy // waves with values
variable p1, p2 // optional range specifiers, default: full wave
If ( !SameNumPnts(wx, wy) )
print "Cov: waves <"+NameOfWave(wx)+"; "+NameOfWave(wy)+"> are not the same length - aborting"
return NAN
elseif ( WaveDims(wx)>1 || WaveDims(wy)>1 )
print "Cov: cannot operate on multidimensional waves - aborting"
endif
p1 = Limit(p1, 0, p1)
p2 = Limit( (ParamIsDefault(p2) ? DimSize(wx,0)-1 : p2), p1, DimSize(wx,0)-1 )
Duplicate/FREE/R=[p1,p2] wx, wxwy
wxwy[0, p2-p1] = wx[ p1+p ] * wy[ p1+p ]
variable foo = (mean(wxwy)-mean(wx,pnt2x(wx,p1),pnt2x(wx,p2))*mean(wy,pnt2x(wy,p1),pnt2x(wy,p2)))
return foo
End
// returns cardinal representation (N, SE, WSW) of <inVal> which is wrapped into 0-360
//
// 2011.09 initial release
Function/S D2Cardinal( inVal )
variable inVal
Make/FREE/T outStr = {"N", "NNE", "NE", "ENE", "E", "ESE", "SE", "SSE", "S", "SSW", "SW", "WSW", "W", "WNW", "NW", "NNW"}
SetScale/P x, 0, 22.5, "", outStr
return outStr( ModWD(inVal) ) // not threadsafe, string
End
// returns <inVal> in radians
//
// 2010.old initial release
ThreadSafe Function D2R( inVal )
variable inVal
return inVal*(PI/180)
End
// returns daqfactory timestamp converted to igor date/time
//
// DAQfactory internally stores time as # of secs since 1970 (01/01?)
//
// 2011.11.07 written
ThreadSafe Function daqfactory2secs( timeval )
variable timeval
return timeval + date2secs(1970, 1, 1)
End
// returns day of the week (1=Sun; ... 7=Sat)
//
// Jan 1, 1904 was a Friday
//
// 2011.11.17 written
Function DayOfWeek( tstamp )
variable tstamp // igor date/time value sec
return mod( 5+trunc(tstamp/86400), 7 )+1
end
// returns density of air in g/m^3 or mol/m^3 using ideal gas law
//
// PV=nRT >> (n/V) = P/(RT) (n/V) = air density mol / m^3 OR g / m^3
// P = barometric presssure mbar
// R = gas constant (mb m^3)/(mol K) OR (mb m^3)/(g K)
// T = temp* Celcius
// *use ambient temp for dry air density, virtual temp for moist air density
//
// 2011.11.15 standarized input units
// 2011.11.09 written
Function DensityOfAir( T_, P_ [, inMoles] )
variable T_ // ambient or virtual temp. Celcius
variable P_ // barometric press. mbar = hPa
variable inMoles // nonzero to return molar density
If ( inMoles )
return P_ / (kRu * (T_+273.15))
else
return P_ / (kR_dryair * (T_+273.15))
endif
end
// applies soft spike filter to entire wave following method of HaPe Schmid; input wave IS modified; optional
// point boundaries permit applying filter to subsection of wave => see IntervalDespikeHaPe() for application to time series
//
// Schmid, HaPe, C. Susan B. Grimmond, Ford Cropley, Brian Offerle, and Hong-Bing Su. "Measurements
// of CO2 and energy fluxes over a mixed hardwood forest in the mid-western United States." Agricultural
// and Forest Meteorology. 103 (2000): 357-374.
// "For each 15min period and variable, the means and variances are calculated. From these diagnostics,
// a threshold for spikes is determined as a multiple of the standard deviation (3.6 S.D. initially, increased
// by 0.3 after each pass). On each pass, a soft spike is registered if the fluctuation from the mean is larger
// than the threshold value, and if the duration of the spike is three or fewer records, corresponding to a
// persistence of 0.3s, for the 10Hz sampling rate. Longer-lasting departures from the period mean are taken
// to indicate possible physical events. After each pass, if spikes are detected, the mean and variance are
// adjusted to exclude data marked as spikes and the process repeated, until either there are no more new
// spikes or the maximum or three iterations is completed (which is rarely the case)."
//
// 2012.06.12 added optional parameter to specify max persistence of registered spikes
// 2011.11.21 split logic to create an Interval_ function and independent function
// 2011.11.18 adapted from legacy function
Function/WAVE DespikeHaPe( wname, tstamp [, multiplier, increment, passes, duration, p1, p2] )
wave wname // target wave ref
wave/Z/D tstamp // double-precision Igor date-time wave; used to evalute spike persistence
// time; if wave does not exist, X-scaling of <wname> is used instead
variable multiplier // optionally specify different base multipier (def: 3.6)
variable increment // optionally specify different per-pass increment (def: 0.3)
variable passes // optionally specify number of passes (def: 3)
variable duration // optionally specify max persistence to be spike, seconds (def: 0.3)
variable p1, p2 // optional inclusive point boundaries
If ( !WaveType(wname) || WaveDims(wname)>1 )
print "DespikeHaPe: source wave <"+NameOfWave(wname)+"> is non-numeric or multidimensional - aborting"
endif
If ( !WaveExists(tstamp) )
Make/D/FREE/N=(DimSize(wname,0)) tstamp = leftx(wname) + deltax(wname)*p
endif
multiplier = ParamIsDefault(multiplier) ? 3.6 : multiplier
increment = ParamIsDefault(increment) ? 0.3 : increment
passes = ParamIsDefault(passes) ? 3 : passes
duration = ParamIsDefault(duration) ? 0.3 : duration
p1 = Limit(p1, 0, p1)
p2 = Limit( (ParamIsDefault(p2) ? DimSize(wname,0)-1 : p2), p1, DimSize(wname,0)-1 )
variable pn, avg, sdev, found, mult, threshold, pp, pp0
Make/FREE/N=(passes,2) results = 0
SetDimLabel 1, 0, spikes, results
SetDimLabel 1, 1, points, results
for (pn=0; pn<passes; pn+=1) // for each pass through time series
avg = mean(wname, pnt2x(wname,p1), pnt2x(wname,p2))
sdev = sqrt( variance(wname, pnt2x(wname,p1), pnt2x(wname,p2)) )
found = 0
mult = multiplier + pn*increment
threshold = abs(mult*sdev)
for (pp=p1; pp<p2+1; pp+=1) // for each point in subinterval
If ( numtype(wname[pp]) ) // for NANs, just skip
continue
elseif ( abs(wname[pp] - avg) > threshold ) // if spike criterea is met
pp0 = found ? pp0 : pp // capture pp=>pp0 only when new spike observed
found = found ? found : 1 // register spike is observed if not already
elseif ( found ) // criterea no longer met but spike was observed
found = 0
If ( (tstamp[pp-1] - tstamp[pp0]) <= duration )
// print/D "spike registered between",pp0,"and",pp-1
wname[pp0, pp-1] = NAN
results[pn][%spikes] += 1
results[pn][%points] += (pp-pp0)
endif
endif
endfor // each point
endfor // each pass
return results
End
// returns dew point based on vapor pressure
//
// Bolton, D., 1980: The Computation of Equivalent Potential Temperature. Monthly Weather
// Review. Vol 108, 1046-1053.
// Eqn 11: T = ( 243.5*ln(es) - 440.8 ) / ( 19.48 - ln(es) ) es = sat. vapor pressure, mbar
//
// 2011.11.08 written
Function DewPoint( e_ )
variable e_ // H2O vapor pressure mbar = hPa
return (243.5 * ln(e_) - 440.8) / (19.48 - ln(e_)) // Celcius
End
// creates boolean waves denoting presence or absence of diagnostics flags
//
// 2011.11.16 renamed from CPC3776diagnostics
// 2011.10.31 initial release
Function/WAVE DiagnoseCPC3776( diagWord, option )
wave diagWord // diagnostic word from CPC 3776
variable option // set bit true to make additional waves (511 = all)
// - - total of boolean flags UINT8 W_cpc3776_flags
// 0 (1) saturator temp UINT8 W_cpc3776_sat_temp
// 1 (2) condensor temp UINT8 W_cpc3776_cond_temp
// 2 (4) optics temp UINT8 W_cpc3776_opt_temp
// 3 (8) inlet flow rate UINT8 W_cpc3776_in_flow
// 4 (16) aerosol flow rate UINT8 W_cpc3776_aero_flow
// 5 (32) laser power UINT8 W_cpc3776_laser_pwr
// 6 (64) liquid level UINT8 W_cpc3776_liq_level
// 7 (128) concentration UINT8 W_cpc3776_conc_flag
// 8 (256) calibration reminder UINT8 W_cpc3776_cal_remind
// ... [all other bits reserved]
variable i, val, n = numpnts(diagWord)
Make/B/U/FREE/N=(n) W_cpc3776_flags, W_cpc3776_sat_temp, W_cpc3776_cond_temp, W_cpc3776_opt_temp
Make/B/U/FREE/N=(n) W_cpc3776_in_flow, W_cpc3776_aero_flow, W_cpc3776_laser_pwr, W_cpc3776_liq_level
Make/B/U/FREE/N=(n) W_cpc3776_conc_flag, W_cpc3776_cal_remind
variable npnts = numpnts(diagWord)
for (i=0; i<npnts; i+=1)
val = diagWord[i]
if ( numtype(val) )
continue
endif
W_cpc3776_sat_temp[i] = TestBit(val, 0) // saturator temp
W_cpc3776_cond_temp[i] = TestBit(val, 1) // condensor temp
W_cpc3776_opt_temp[i] = TestBit(val, 2) // optics temp
W_cpc3776_in_flow[i] = TestBit(val, 3) // inlet flow rate
W_cpc3776_aero_flow[i] = TestBit(val, 4) // aerosol flow rate
W_cpc3776_laser_pwr[i] = TestBit(val, 5) // laser power
W_cpc3776_liq_level[i] = TestBit(val, 6) // liquid level
W_cpc3776_conc_flag[i] = TestBit(val, 7) // concentration flag
W_cpc3776_cal_remind[i] = TestBit(val, 8) // calibration reminder
W_cpc3776_flags[i] = W_cpc3776_sat_temp[i]+W_cpc3776_cond_temp[i]+W_cpc3776_opt_temp[i]+W_cpc3776_in_flow[i]
W_cpc3776_flags[i] += W_cpc3776_aero_flow[i]+W_cpc3776_laser_pwr[i]+W_cpc3776_liq_level[i]
W_cpc3776_flags[i] += W_cpc3776_conc_flag[i]+W_cpc3776_cal_remind[i]
endfor
Duplicate/B/U/O W_cpc3776_flags, :W_cpc3776_flags
If ( TestBit(option, 0) )
Duplicate/B/U/O W_cpc3776_sat_temp, :W_cpc3776_sat_temp
endif
If ( TestBit(option, 1) )
Duplicate/B/U/O W_cpc3776_cond_temp, :W_cpc3776_cond_temp
endif
If ( TestBit(option, 2) )
Duplicate/B/U/O W_cpc3776_opt_temp, :W_cpc3776_opt_temp
endif
If ( TestBit(option, 3) )
Duplicate/B/U/O W_cpc3776_in_flow, :W_cpc3776_in_flow
endif
If ( TestBit(option, 4) )
Duplicate/B/U/O W_cpc3776_aero_flow, :W_cpc3776_aero_flow
endif
If ( TestBit(option, 5) )
Duplicate/B/U/O W_cpc3776_laser_pwr, :W_cpc3776_laser_pwr
endif
If ( TestBit(option, 6) )
Duplicate/B/U/O W_cpc3776_liq_level, :W_cpc3776_liq_level
endif
If ( TestBit(option, 7) )
Duplicate/B/U/O W_cpc3776_conc_flag, :W_cpc3776_conc_flag
endif
If ( TestBit(option, 8) )
Duplicate/B/U/O W_cpc3776_cal_remind, :W_cpc3776_cal_remind
endif
wave outref = W_cpc3776_flags
return outref
End
// creates boolean waves denoting presence or absence of diagnostics flags
//
// 2013.05.23 made arg `option` optional
// 2011.11.16 renamed from CSAT3diagnostics()
// 2011.10.31 reduced If(...) structure to boolean additions under default switch statment
// 2011.09 initial release
Function/WAVE DiagnoseCSAT3( diagWord [, option] )
wave diagWord // diagnostic word from CSAT3
variable option // set bit true to make additional waves (15=all)
// 0 (1) speed of sound diff. boolean W_csat_del_T
// 1 (2) poor signal lock boolean W_csat_sig_lck
// 2 (4) high amplitude boolean W_csat_amp_h
// 3 (8) low amplitude boolean W_csat_amp_l
// ... [all other bits reserved]
variable i, val, n = numpnts(diagWord)
option = (ParamIsDefault(option) ? 0 : option)
Make/B/U/FREE/N=(n) W_csat3_flags, W_csat3_del_T, W_csat3_sig_lck, W_csat3_amp_h, W_csat3_amp_l
for (i=0; i<numpnts(diagWord); i+=1)
val = diagWord[i]
if ( numtype(val) )
continue
endif
switch (val)
case -99999:
case 61440:
case 61503:
case 61441:
case 61442:
W_csat3_flags[i] += 1
break;
default:
W_csat3_del_T[i] = TestBit(val, 15) // speed of sound diff
W_csat3_sig_lck[i] = TestBit(val, 14) // poor signal lock
W_csat3_amp_h[i] = TestBit(val, 13) // high amp
W_csat3_amp_l[i] = TestBit(val, 12) // low amp
W_csat3_flags[i] = W_csat3_del_T[i]+W_csat3_sig_lck[i]+W_csat3_amp_h[i]+W_csat3_amp_l[i]
endswitch
endfor
Duplicate/B/U/O W_csat3_flags, :W_csat3_flags
If ( TestBit(option, 0) )
Duplicate/B/U/O W_csat3_del_T, :W_csat3_del_T
endif
If ( TestBit(option, 1) )
Duplicate/B/U/O W_csat3_sig_lck, :W_csat3_sig_lck
endif
If ( TestBit(option, 2) )
Duplicate/B/U/O W_csat3_amp_h, :W_csat3_amp_h
endif
If ( TestBit(option, 3) )
Duplicate/B/U/O W_csat3_amp_l, :W_csat3_amp_l
endif
wave outref = W_csat3_flags
return outref
End
// creates boolean waves denoting presence or absence of diagnostics flags
//
// 2011.11.16 renamed from LI7500diagnostics()
// 2011.10.31 initial release
Function/WAVE DiagnoseLI7500( diagWord, option )
wave diagWord // diagnostic word from LI-7500
variable option // set bit true to make additional waves (31=all)
// - - total of boolean flags UINT8 W_li7500_flags
// 0 (1) chopper warning UINT8 W_li7500_chopper
// 1 (2) detector warning UINT8 W_li7500_detector
// 2 (4) phase lock loop warning UINT8 W_li7500_pll
// 3 (8) sync warning UINT8 W_li7500_sync
// 4 (16) AGC value SP W_li7500_AGC
// ... [all other bits reserved]
variable i, val, n = numpnts(diagWord)
Make/B/U/FREE/N=(n) W_li7500_flags, W_li7500_chopper, W_li7500_detector, W_li7500_pll, W_li7500_sync
Make/FREE/N=(n) W_li7500_AGC
for (i=0; i<numpnts(diagWord); i+=1)
val = (diagWord[i] %^ 0xF0 ) // swap bits 4-7: 1111 0000
if ( numtype(val) )
continue
endif
W_li7500_chopper[i] = TestBit(val, 7) // chopper not OK
W_li7500_detector[i] = TestBit(val, 6) // detector not OK
W_li7500_pll[i] = TestBit(val, 5) // phase lock loop not OK
W_li7500_sync[i] = TestBit(val, 4) // sync not OK
W_li7500_flags[i] = W_li7500_chopper[i]+W_li7500_detector[i]+W_li7500_pll[i]+W_li7500_sync[i]
W_li7500_AGC[i] = ( val & 0x0F ) * 6.25
endfor
Duplicate/B/U/O W_li7500_flags, :W_li7500_flags
If ( TestBit(option, 0) )
Duplicate/B/U/O W_li7500_chopper, :W_li7500_chopper
endif
If ( TestBit(option, 1) )
Duplicate/B/U/O W_li7500_detector, :W_li7500_detector
endif
If ( TestBit(option, 2) )
Duplicate/B/U/O W_li7500_pll, :W_li7500_pll
endif
If ( TestBit(option, 3) )
Duplicate/B/U/O W_li7500_sync, :W_li7500_sync
endif
If ( TestBit(option, 4) )
Duplicate/O W_li7500_AGC, :W_li7500_AGC
endif
wave outref = W_li7500_flags
return outref
End
// averages values by day; does not handle NANs yet
//
// 2011.11.17 added <mode> to omit days if desired
// 2011.11.09 written
Function/WAVE DielAverage( wname, tstamp, mode)
wave wname // value wave
wave/D tstamp // time wave
variable mode // bitmask: set bit HI to exclude that day of the week
// weekends only = 2+4+8+16+32 = 62
// weekdays only = 1+64 = 65
// bit value DOW omitted bit value DOW omitted
// - 0 -none- 3 8 Wednesday
// 0 1 Sunday 4 16 Thursday
// 1 2 Monday 5 32 Friday
// 2 4 Tuesday 6 64 Saturday
variable i, lo, hi, oi
If ( mode )
for (i=0; i<7; i+=1)
If ( TestBit( mode, i ) )
Extract/FREE wname, tmp0, DayOfWeek(tstamp) != (i+1)
Extract/D/FREE tstamp, tmp1, DayOfWeek(tstamp) != (i+1)
wave wname = tmp0
wave/D tstamp = tmp1
endif
endfor
endif
Duplicate/D/FREE tstamp, tdiel
tdiel = mod(tstamp, 86400)
Duplicate/FREE wname, sortw
Sort tdiel, tdiel, sortw
Make/FREE/N=( 86400 / (tstamp[1]-tstamp[0]) ) dielw = NAN
SetScale/P x, 0, (tstamp[1]-tstamp[0]), "dat", dielw
for (i=0; 1; i+=1)
if ( i+1 >= DimSize(tdiel,0) )
dielw[oi] = mean(sortw, pnt2x(sortw,lo), pnt2x(sortw, i ) )
break
elseif ( tdiel[i+1] > tdiel[i] )
hi = i
dielw[oi] = mean(sortw, pnt2x(sortw,lo), pnt2x(sortw,hi) )
lo = i+1
oi += 1
endif
endfor
return dielw
End
// converts decimal day of year into igor date/time
Function doy2sec( doy, year )
variable doy // day of year, indexed from 1
variable year // year is reqd since doy depends on year
year = trunc(year)
Make/FREE/N=366 monthmap = NAN; // monthmap[0] not used
monthmap[1,31]=1; monthmap[32,59]=2; monthmap[60,90]=3
monthmap[91,120]=4; monthmap[121,151]=5; monthmap[152,181]=6
monthmap[182,212]=7; monthmap[213,243]=8; monthmap[244,273]=9
monthmap[274,304]=10; monthmap[305,334]=11; monthmap[335,365]=12
Make/FREE/N=13 endDOY = {NAN, 31, 59, 90, 120, 151, 181, 212, 242, 273, 304, 334, 365}
If ( !mod(year,4) && ( mod(year,100) || !mod(year,400) ) ) // LEAP YEAR ADJ
// if (year % 4 == 0) AND ((year % 100 !=0) OR (year % 400 == 0))
InsertPoints 60, 1, monthmap;
monthmap[60] = 2
endDOY[2, 12] += 1
doy = mod(doy, 366)
else
doy = mod(doy, 365)
endif
variable mo = monthmap[trunc(doy)]
return date2secs(year, mo, (doy - endDOY[mo-1])) + mod(doy,1)*86400
End
// TODO
// verify this function
//
// returns latent heat flux using eddy covariance
//
// 2011.11.11 written
Function ECLatentHeat( h2o, w_ [, T_, p1, p2] )
wave h2o // water vapor molar density g / m^3
wave w_ // vertical wind m/s
wave T_ // optional ambient temp Celcius
// temp determines latent heat of vaporization as f(meanT); uses 2.5e3 J/g otherwise
variable p1, p2 // optional point boundaries, inclusive
If ( !SameNumPnts(h2o, w_) || (WaveExists(T_) && !SameNumPnts(h2o, T_)) )
print "ECLatentHeat: input waves were not the same length - aborted"
return NAN
endif
p1 = Limit(p1, 0, p1)
p2 = Limit( (ParamIsDefault(p2) ? DimSize(h2o,0) : p2), p1, DimSize(h2o,0) )
variable out = Cov( h2o, w_, p1=p1, p2=p2 )
If ( ParamIsDefault(T_) || !WaveExists(T_) )
out *= 2.5e3
else
out *= LatentHeatVapH2O(mean(T_,pnt2x(T_,p1),pnt2x(T_,p2)))
endif
return out
End
// returns sensible heat flux in Watts/meter^2 using eddy covariance
//
// weak source: Eqn 11.20 Arya, S. Pal. Introduction to Micrometerology. 2nd Ed. 2001. Academic Press.
//
// 2013.05.20 fix numpnts check
// 2012.06.18 *bug* fixed units scaling which previously resulted in values 100X too small
// 2011.11.11 written
Function ECSensibleHeat( T_, w_ , P_, Q_, [p1, p2] )
wave T_ // ambient temp. Celcius
wave w_ // vertical wind m/s
wave P_ // ambient press. mbar
wave Q_ // specific humidity dimensionless, 0-1
variable p1, p2 // optional point boundaries, inclusive
Make/FREE/WAVE/N=5 wlist = {T_, w_, P_, Q_}
If ( !SameNumPntsW( wlist ) )
print "ECSensibleHeat: input waves were not the same length - aborted"
return NAN
endif
p1 = Limit(p1, 0, p1)
p2 = Limit( (ParamIsDefault(p2) ? DimSize(T_,0) : p2), p1, DimSize(T_,0) )
variable out, mt, mq, mp
mt = mean(T_, pnt2x(T_,p1), pnt2x(T_,p2) )
mq = mean(Q_, pnt2x(Q_,p1),pnt2x(Q_,p2) )
mp = mean(P_, pnt2x(P_,p1),pnt2x(P_,p2) )
out = 100 * kCp_dryair * (1+0.84*mq) * DensityOfAir( VirtualTemp( mt, mq ), mp ) * Cov( T_, w_, p1=p1, p2=p2 )
return out
End
// return eddy covariance flux of CO2 from gas density and wind measurements
// value will be NAN if arguments contain NAN
//
// 2013.05.20 written
Function EC_co2( co2, w_, [p1, p2] )
wave co2 // density of carbon dioxide g / m^2 or mol/m^2
wave w_ // vertical wind speed m / s
variable p1, p2 // optional point boundaries, inclusive
If ( !SameNumPnts(co2, w_) )
print "EC_co2: input waves are not all the same length - aborting"
return NAN
endif
p1 = Limit(p1, 0, p1)
p2 = Limit( (ParamIsDefault(p2) ? DimSize(co2,0)-1 : p2), p1, DimSize(co2,0)-1 )
return Cov( co2, w_, p1=p1, p2=p2 )
End
// TODO
// verify results of this function
//
// returns eddy covariance density corrections following WPL(80) using variables
// value returned is complex: water vapor term resides in real portion, temp. term in imaginary portion
//
// <Fc> = <Fco> + <RHOc> * [ mu*Cov( W, RHOv ) / <RHOd> + (1+ mu*sigma)*Cov( W, T ) / <T> ]
// |_mult_| |___water vapor term______| |___temperature term_________|
//
// Fco uncorrected flux g m-2 s-1
// RHOc density of constituent C g m-3
// mu (MW dry air / MW water vapor) = 1/0.622 dimensionless
// W vertical wind speed m s-1
// RHOv density of water vapor g m-3
// RHOd density of dry air g m-3
// sigma = RHOv / RHOd dimensionless
// T temperature Kelvin
// <...> denotes averaging operation
//
// 2011.11.22 written
Function/C EC_WPL80( meanRHOc, meanRHOv, meanRHOd, meanT, cov_w_rhoV, cov_w_T )
variable meanRHOc // mean mass density of constituent C g / m^3
variable meanRHOv // mean mass density of water vapor g / m^3
variable meanRHOd // mean mass density of dry air g / m^3
variable meanT // mean temperature Celcius
variable cov_w_rhoV // covariance of vertical wind and rhoV g m-2 s-1
variable cov_w_T // covariance of vertical wind and temp C m / s = K m / s
variable mu = kMW_dryair / kMW_h2o
variable W_corr = mu * cov_w_rhoV / meanRHOd
variable T_corr = (1 + mu*(meanRHOv / meanRHOd)) * cov_w_T / meanT
return cmplx( meanRHOc*W_corr, meanRHOc*T_corr )
End
// TODO
// verify results of this function
//
// returns eddy covariance density corrections following WPL (1980) using time series
// value returned is complex: water vapor term resides in real portion, temp. term in imaginary portion
//
// <Fc> = <Fco> + <RHOc> * [ mu*Cov( W, RHOv ) / <RHOd> + (1+ mu*sigma)*Cov( W, T ) / <T> ]
// |_mult_| |___water vapor term______| |___temperature term_________|
//
// Fco uncorrected flux g m-2 s-1
// RHOc density of constituent C g m-3
// mu (MW dry air / MW water vapor) = 1/0.622 dimensionless
// W vertical wind speed m s-1
// RHOv density of water vapor g m-3
// RHOd density of dry air g m-3
// sigma = RHOv / RHOd dimensionless