IDL Script RFH_NCURVE.idl
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; FILE : /fs/cgd/home0/davestep/idl_routines/RFH_NCURVE.idl
; AUTHOR : David Stepaniak, NCAR/CGD/CAS
; DATE INITIATED : 15 September 1998
; LAST MODIFIED : Wed Sep 16 12:28:58 MDT 1998
;
; DESCRIPTION : Generates a relative frequency histogram for a given dis-
; crete (or nearly continuous) and presumably random variable
; (perhaps such as surface temperature anomalies), and overlays
; a normal curve with the same standard deviation as the vari-
; able and same area as its histogram. The data of the input
; variable is transformed to have a mean of 0, and the x-axis
; is plotted in units of standard deviation.
;
; NOTE : Script may be run by invoking IDL from the command line:
; e.g. /fs/cgd/home0/davestep/idl_routines/idl
; Then typing '.r RFH_NCURVE.idl' at the IDL prompt:
; IDL>.r RFH_NCURVE.idl
;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; REQUIRED INPUT PARAMETERS SUPPLIED BY USER:
file = 'Boulder_mm_ts_1897_1998.32ieee' ; Data file.
; Data is assumed to be written
; as one unformatted sequential
; access Fortran record, e.g.
; WRITE (9) x
; where x is a one-dimensional
; array.
ndata = 102*12 ; Number of data points, including
; those data points assigned a mis-
; sing vlaue.
msg = -99.9 ; Missing value. If there are no
; data with the missing value, then
; any number of type float will do.
; The following three parameters
; must have same units as the data.
; Also, data will be transformed
; so that the mean is zero; thus,
; symmetric left and right center
; bin limits are recommended.
ctd_lftm_bin = -10. ; Center point of the leftmost bin.
ctd_rgtm_bin = 10. ; Center point the rightmost bin.
dbin = 0.2 ; Binsize. (Spacing between center
; points of adjacent bins.)
; (The range between ctd_lftm_bin
; and ctd_rgtm_bin should be di-
; visible by dbin without
; remainder.
; Plot annotations:
; ----------------
; Plot title:
p_title = 'Boulder Monthly Surface Temperature Anomalies!C' + $
'1897-1998 (1961-90 Base Period)'
; Format of standard deviation to be displayed on plot,
; f77/f90 formatting conventions:
p_sd_format = '(F4.2)'
; Units of data:
data_units = '!Eo!NC'
; Miscellaneous:
; --------------
; Number of minor tickmarks per
; major tickmark interval:
p_nxminor = 4
; Plot axes thickness and tickmark thickness factor (default
; is 1.0):
p_thick = 2.5
; Font size factor (default is 1.0):
p_fsize = 1.1
; Font thickness factor(default is 1.0):
p_fthick = 2.5
; Histogram and normal curve line thickness factor (default is
; 1.0):
p_lthick = 2.5
; PostScript options:
; -------------------
; PostScript file output option, 'y'
; for affirmative, 'n' for default X
; display (if Postcript output, file
; will be in portrait mode with name
; 'RFH_NCURVE.ps'.)
PS = 'y'
; If PS = 'n' then following two options may be ignored:
; Width of PS plot in inches (including annotations):
PSw = 5.0
; Height of PS plot in inches (including annotations):
PSh = 5.0
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Data file IO:
data = FLTARR(ndata)
OPENR, 9, file, /F77_UNFORMATTED
READU, 9, data
CLOSE, 9
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Account for missing value, put all nonmissing data in the one-dimensional
; array x of size n, compute mean and standard deviation of x:
n = 0
FOR i = 0, ndata - 1 DO BEGIN
IF ( data(i) ne msg ) THEN BEGIN
n = n + 1
ENDIF
ENDFOR
x = FLTARR(n)
k = -1
FOR i = 0, ndata - 1 DO BEGIN
IF ( data(i) ne msg ) THEN BEGIN
k = k + 1
x(k) = data(i)
ENDIF
ENDFOR
mn_var_01 = MOMENT(x)
mn = mn_var_01(0)
sd = SQRT( mn_var_01(1) )
; Convert numeric value of standard deviation to string for plotting
; purposes:
p_sd_string = STRING(FORMAT=p_sd_format,sd)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Transform x so that its mean is 0 (DO NOT normalize by standard deviation):
x = x - mn
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Compute area of RELATIVE frequency histogram of x using bin
; input information. (The area should equal dbin.)
n_bin = ROUND( (ctd_rgtm_bin - ctd_lftm_bin) / dbin ) + 1
; Number of bins.
n_bin_lims = n_bin + 1
; Number of bin boundaries.
bin_lims = (ctd_lftm_bin - 0.5*dbin) + FINDGEN(n_bin_lims)*dbin
; Positions of bin boundaries.
c_bins = ctd_lftm_bin + FINDGEN(n_bin)*dbin
; Positions of bin centers.
ARFH = 0.
; Initialize area of relative frequency histogram as 0.
FOR i = 1, n_bin_lims - 1 DO BEGIN
p = 0
FOR k = 0, n - 1 DO BEGIN
IF ( ( x(k) gt bin_lims(i-1) ) and ( x(k) le bin_lims(i) ) ) THEN BEGIN
p = p + 1
ENDIF
ENDFOR
ARFH = ARFH + ( FLOAT(p) / FLOAT(n) )*dbin
ENDFOR
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Compute percentage of data which lie between (+/-)1 sd:
p = 0
FOR k = 0, n - 1 DO BEGIN
IF ( ( x(k) ge -sd ) and ( x(k) le sd ) ) THEN BEGIN
p = p + 1
ENDIF
ENDFOR
pdpm1sd = FLOAT(p)*100./FLOAT(n)
; Convert numeric value of percentage to string for plotting purposes:
pdpm1sd_string = STRING(FORMAT='(F4.1)',pdpm1sd)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Generate relative frequency histogram:
H = HISTOGRAM(x, BINSIZE=dbin, MIN=ctd_lftm_bin, MAX=ctd_rgtm_bin)/FLOAT(n)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Generate normal curve with same sd as x and same area as the
; relative frequency histogram:
y = ( ARFH / ( sd*SQRT(2.*!PI) ) ) * EXP( - ((c_bins*c_bins)/(sd*sd))/2 )
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Compute x-axis tickmark parameters so that x-axis is normalized to
; standard deviation:
; Determine number of major tickmark intervals:
p_xticks_float = (ctd_rgtm_bin - ctd_lftm_bin)/sd
p_xticks = CEIL(p_xticks_float)
IF ( ( p_xticks MOD 2 ) ne 0 ) THEN BEGIN
p_xticks = p_xticks + 1
ENDIF
; Compute major tickmark values:
n_xtickv = p_xticks + 1
p_xtickv = (-FLOAT(p_xticks/2) + FINDGEN(n_xtickv)*1.)*sd
; Assign major tickmark labels to reflect normalization of x-axis to sd:
p_xtickname = STRARR(n_xtickv)
FOR k = 0, n_xtickv - 1 DO BEGIN
tm_int_form = -FLOAT(p_xticks/2) + k
IF ( tm_int_form lt 0 ) THEN BEGIN
p_xtickname(k) = STRING(FORMAT='(I2)',tm_int_form)
ENDIF
IF ( tm_int_form ge 0 ) THEN BEGIN
p_xtickname(k) = STRING(FORMAT='(I1)',tm_int_form)
ENDIF
ENDFOR
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; PostScript or X display:
IF (PS eq 'y') THEN BEGIN
xsz = PSw
ysz = PSh
xof=(8.5-xsz)/2.
yof=(11.0-ysz)/2.
SET_PLOT, 'PS'
DEVICE, /PORTRAIT, FILE='RFH_NCURVE.ps', /INCHES, XOFFSET=xof, YOFFSET=yof, $
XSIZE=xsz, YSIZE=ysz
ENDIF
IF (PS eq 'n') THEN BEGIN
WINDOW, XSIZE = 600, YSIZE = 700, TITLE = 'DATA HISTOGRAM AND NORMAL CURVE'
ENDIF
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Plot histogram and normal curve:
PLOT, c_bins, H, PSYM = 10, $
TITLE = p_title, $
XTITLE = 'Standard Deviations (1 sd = ' + $
p_sd_string + data_units + ')', $
SUBTITLE = '(' + pdpm1sd_string + '% of data between !9+!31 sd)', $
YTITLE = 'Relative Frequency', $
XTICKS = p_xticks, $
XTICKV = p_xtickv, $
XTICKNAME = p_xtickname, $
XMINOR = p_nxminor + 1, $
XTHICK = p_thick, $
YTHICK = p_thick, $
CHARSIZE = p_fsize, $
CHARTHICK = p_fthick, $
THICK = p_lthick
OPLOT, c_bins, y, THICK = p_lthick
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Do actual PS output (PS file must be closed to do so), reset to X:
IF (PS eq 'y') THEN BEGIN
DEVICE,/PORTRAIT,/CLOSE
SET_PLOT, 'X'
ENDIF
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
END
Back