Changeset 108

Timestamp:

06/03/08 03:15:28 (16 years ago)

Author:

ericg

Message:

Bug in computation of 1D stddev in trends.pro corrected

File:

• trunk/procs/trends.pro (modified) (5 diffs)

trunk/procs/trends.pro

@@ -127 +127 @@
                ENDIF 
                IF fld_flag EQ 2 THEN fldrem_t2 = fldrem
+               sc_ampl = sqrt((moment(fldrem))[1])
+               sc_ampl2 = max(fldrem)-min(fldrem)
                fldrem = reform(fldrem#replicate(1, long(lenght/running)), lenght)
 
@@ -133 +135 @@
                err_std = sqrt( total( (temperr-mean(fld))^2 )/n_elements(temperr) )
                
-               IF NOT keyword_set(silent) THEN print, '     1D times-serie info         : mean +- err ', (moment(fld))[0], err_std
+               IF NOT keyword_set(silent) THEN BEGIN
+                  print, '     1D times-serie info         : mean +- err ', (moment(fld))[0], err_std
+                  print, '       Mean seasonal cycle amplitude stdev / max-min', sc_ampl,  sc_ampl2
+               ENDIF 
 
                ; Compute std error-bar using running bootstrap window (width in plt_def)
                IF jpt GT boot_win(0) THEN BEGIN 
                   temperr = stat_error(fld, boot_win(0))
-                  err_std_1 = sqrt( total((temperr-sqrt((moment(fld))[1]) )^2)/n_elements(temperr) )
+                  err_std_1 = sqrt( total((temperr-sqrt((moment(fld - fldrem))[1]) )^2)/n_elements(temperr) )
                ENDIF ELSE err_std_1 = 0
                IF jpt GT boot_win(1) THEN BEGIN 
                   temperr = stat_error(fld, boot_win(1))
-                  err_std_2 = sqrt( total((temperr-sqrt((moment(fld))[1]) )^2)/n_elements(temperr) )
+                  err_std_2 = sqrt( total((temperr-sqrt((moment(fld - fldrem))[1]) )^2)/n_elements(temperr) )
                ENDIF  ELSE err_std_2 = 0
                IF jpt GT boot_win(2) THEN BEGIN 
                   temperr = stat_error(fld, boot_win(2))
-                  err_std_3 = sqrt( total((temperr-sqrt((moment(fld))[1]) )^2)/n_elements(temperr) )
+                  err_std_3 = sqrt( total((temperr-sqrt((moment(fld - fldrem))[1]) )^2)/n_elements(temperr) )
                ENDIF  ELSE err_std_3 = 0
                
-               IF NOT keyword_set(silent) THEN print, '     1D times-serie anomaly info : var/stddev +- error (24,36,48)', (moment(fld))[1], sqrt((moment(fld))[1]), err_std_1, err_std_2, err_std_3 
+               IF NOT keyword_set(silent) THEN print, '     1D times-serie anomaly info : var/stddev +- error (24,36,48)', (moment(fld - fldrem))[1], sqrt((moment(fld - fldrem))[1]), err_std_1, err_std_2, err_std_3 
                IF NOT keyword_set(silent) THEN print, ' '
             ENDIF ELSE fldrem =  mean(fld)
@@ -236 +241 @@
             ENDCASE
 
-;            print, min(fld[where(fld le valmask/10)]), max(fld[where(fld le valmask/10)])
-;            print, min(fldrem[where(fldrem le valmask/10)]), max(fldrem[where(fldrem le valmask/10)])
              idx = where(fld EQ valmask)
              fld = fld - fldrem
@@ -243 +246 @@
              mean_sc = fldrem
              IF idx(0) NE -1 THEN fld(idx) = valmask
-;            print, min(fld[where(fld le valmask/10)]), max(fld[where(fld le valmask/10)])
          ENDELSE  
 
@@ -369 +371 @@
          flddesc.data = fld
          cdf_description = nc_build(cmd_wrk, flddesc, type, vargrid) 
-         ;;fldcdf = {data:fld, units:field.units, short_name:cmd_wrk.var, long_name:field.legend+' ('+legbox+')', missing_value:valmask, direc:type}
          fldcdf = {data:fld, units:field.units, short_name:field.legend+'_'+legbox, long_name:field.legend+' ('+legbox+')', missing_value:valmask, direc:type}
          file_out = cmd_wrk.var+'_'+strtrim(string(FORMAT = '(I2.2)', ABS(write_data)), 2)+'.nc'