;*********************************************************************************************
;Spectrum/spectra calculation and plotting using the precipitation anomaly
;*********************************************************************************************

begin

;*********************************************************************************************
;Monthly precipitation anomaly(1901-2018)
;*********************************************************************************************
;f    = addfile("stdanom2_1.nc","r")
;appt = f->pre(:,0,0)

appt= asciiread("cru_wavelet_det.txt",-1,"float")
;appt@_FillValue = getFillValue(appt)

;*********************************************************************************************
;detrending opt: 0=>remove the mean and 1=>remove mean and trend
;*********************************************************************************************
d = 0

;*********************************************************************************************
;Smoothing of the periodogram. Must be odd and minimum value of 3
;*********************************************************************************************
sm = 3

;*********************************************************************************************
;Smoothing percentage: (0.0<=pct<=1.0) with a value of 0.10 being the most common
;*********************************************************************************************
pct = 0.10

;*********************************************************************************************
;calculate spectrum
;*********************************************************************************************
sdof = specx_anal(appt,d,sm,pct)

printVarSummary(sdof)

splt = specx_ci(sdof,0.05,0.95)                ;calculates confidence interval of the red noise

;Degrees of freedom returned from specx_anal function.
;Lower Confidence Interval Limit (0.0 < lowval <1.).Typical value is 0.05.
;Upper Confidence Interval Limit (0.0 < highval <1.). Typical value is 0.95.

printVarSummary(splt)
;exit
;splt results in 4 curves. These curves are used to generate the graph. Note in the figure that there are 4 curves in the order below:

;splt(0,:) - input spectrum
;splt(1,:) - Markov "Red Noise" spectrum
;splt(2,:) - lower confidence bound of Markov
;splt(3,:) - upper confidence bound of Markov

;*********************************************************************************************
;plotting
;*********************************************************************************************
wks  = gsn_open_wks("png","spec_normal")                                    ; send graphics to PS file 

res               = True                                            ; Enables graphics customization
res@gsnLeftString  = "PSD"
res@tiXAxisString = "Frequency (cycles/year)"
res@tiYAxisString = "Variance"
res@tmXTOn        = False                                            ; Disable top x-axis (XT) minortick.
res@xyLineThicknesses      = (/3.,2.,2.,2./)                         ; Line thicknesses.
res@xyDashPatterns         = (/0,0,1,1/)                             ; Line types.
res@xyLineColors           = (/"foreground","green","blue","red"/)
res@pmLegendDisplayMode    = "Always"                                ; Enables legend.
res@pmLegendSide           = "Top"                                   ; Legend posistion.
res@pmLegendParallelPosF   = 0.2                                    ; Move legend in x direction.
res@pmLegendOrthogonalPosF = -0.5                                  ; Move legend in y direction.
res@pmLegendWidthF         = 0.1                                     ; Legend width.
res@pmLegendHeightF        = 0.15                                    ; Legend height.
res@lgPerimOn              = False                                   ; Desable box around legend.
res@lgLabelFontHeightF     = 0.024                                   ; Legend label font height.
res@lgItemOrder            = (/3,2,1,0/)                             ; Reorder subtitle labels.
res@vpHeightF = 0.43                                                 ; Changes the aspect ratio
 res@vpWidthF  = 1.00

 ; font size
  res@tiMainFontHeightF = 0.022     ; fontsize of the main title
  res@txFontHeightF = 0.02          ; fontsize of the subtitles
  res@tmXBLabelFontHeightF = 0.018  ; fontsize of tickmark labels (x-axis)
  res@tmYLLabelFontHeightF = 0.018  ; fontsize of tickmark labels (x-axis)


res@gsnMaximize             =True
;res@trYLog                  =True                                  ;log scaling
;res@trYMaxF                = 700                                    ; Maximum value of y.
;res@trYMinF                = 0.1                                  ; Minimum value of y.
res@trXMaxF                = 0.5                                   ; Maximum value of x.
res@trXMinF                = 0.0                                   ; Minimum value of x.
;res@tmYLMode               = "Manual"                             ; Y axis customization my way.
;res@tmYLTickStartF         = res@trYMinF
;res@tmYLTickSpacingF       = 2
;res@tmYLTickEndF           = res@trYMaxF
res@tmXBMode               = "Manual"                                ; X  axis customization my way.
res@tmXBTickStartF         = res@trXMinF
res@tmXBTickSpacingF       = 0.1
res@tmXBTickEndF           = res@trXMaxF
res@tmYLMinorOn            = True                                   ; Disable minor tickmark on left y-axis (YL).
res@tmXBMinorOn            = True                                   ; Disable minor tickmark on bottom x-axis (XB).
res@tmYLFormat             = "0@;*.2f"                               ; Formats the left axis (YL) numbers. 
res@xyExplicitLegendLabels = (/"Spectra","Red Noise","5% significance","95% significance"/)                      
                               
plot=gsn_csm_xy(wks,sdof@frq,splt,res)                               ; Generate 4 graph with 4 curves. 

end