;---------------------------------------------------------------------- ; demod_cmplx_1.ncl ; ; Concepts illustrated: ; - Reading a simple text file (Wolf Sunspot Numbers) via 'asciiread' ; - Creating variables for input to 'demod_cmplx' ; - Use 'demod_cmplx' ; - Explicitly extracting the amplitudes and ; phase variables from the retrned variable of type 'list' ; - Use 'unwrap_phase' to make better phase plots ; - Use 'gsn_attach_plots' ;---------------------------------------------------------------------- ; Bloomfield, P. (1976) ; Fourier Analysis of Time series: An Introduction ; Wiley , 1976: Chapter 6 ;-- ; Bloomfield, P. (2000) ; Fourier Analysis of Time series: An Introduction ; Wiley , 2000: Chapter 7 ;---------------------------------------------------------------------- ; Sunspot numbers are Version 1 to match those used by Bloomfield. ; Source: WDC-SILSO, Royal Observatory of Belgium, Brussels ; http://www.sidc.be/silso/datafiles ;---------------------------------------------------------------------- ; All necessary NCL libraries are automatically loaded. ;---------------------------------------------------------------------- diri = "./" fili = "sunspot_wolf.1700_2014.yearssn.txt" ; 1st row (line) has been removed nyrs = numAsciiRow(diri+fili) SSN = asciiread(diri+fili, (/nyrs,2/), "float") SSN@long_name = "Sunspot Number: Source: WDC-SILSO, Royal Observatory of Belgium, Brussels" printVarSummary(SSN) printMinMax(SSN(:,1), 0) print(" ") ; for clarity & convenience, explicitly extract the values yyyy = SSN(:,0) ssn = SSN(:,1) ssn@long_name = "Sun Spot Number" ntim = dimsizes(yyyy) ; # of years ;---Detrend (optional) ;ssn = dtrend(ssn, True) ;ssn@long_name = "Sun Spot Number: detrended" ;---Specify demodulation period/frequency (could be fractional) period = 11.0 ; Bloomfield uses this frqdem = 1.0/period ;---Perform complex demodulation on the anomaly series nwt = 41 ; same # of pts as Bloomfield BUT different filter frqcut = 0.5*frqdem ; arbitrary low-pass cutoff frequency ;;frqcut = 0.025 ; could be explicitly set; try different values ssndm = demod_cmplx(ssn, frqdem, frqcut, nwt, 0, False) print(ssndm) ; type list ;---Convenience and clarity ;---Explicitly extract variable(s) from returned list variable ssnAmp = ssndm[0] ; [0] list syntax; Amp(time) ssnPha = ssndm[1] ; [1] Pha(time) ssnPhau = ssndm[2] ; [2] Pha_Unwrapped(time) delete(ssndm) ; no longer needed printVarSummary(ssnAmp) printMinMax(ssnAmp,0) print(" ") printVarSummary(ssnPha) printMinMax(ssnPha,0) print(" ") printVarSummary(ssnPhau) printMinMax(ssnPhau,0) print(" ") ;=============================================================== ; PLOT ;====================================== yrStrt = toint(yyyy(0)) yrLast = toint(yyyy(ntim-1)) plot = new (4, graphic) wks = gsn_open_wks ("png","demod_cmplx") ; send graphics to PNG file res = True ; plot mods desired res@gsnDraw = False ; don't draw frame yet res@gsnFrame = False ; don't advance frame yet res@vpHeightF= 0.4 ; change aspect ratio of plot res@vpWidthF = 0.85 res@vpXF = 0.100 ; move left edge res@trXMinF = yrStrt res@trXMaxF = yrLast+1 res@tiYAxisString = "SSN" ; ssn@long_name ; y-axis label res@tiMainString = "Sun Spot Number:"+yrStrt+"-"+yrLast plot(0) = gsn_csm_xy (wks,yyyy,ssn,res) delete(res@tiYAxisString) plot(1) = gsn_csm_xy (wks,yyyy,ssnAmp ,res) getvalues plot(1) ; add a information label to amplitude plot "tmYLLabelFontHeightF" : fheight end getvalues txres = True txres@txFontHeightF = fheight ; Match font height label = "p="+sprintf("%3.1f", period)+": fdem="+sprintf("%5.3f", frqdem) \ +": fcut="+sprintf("%5.3f", frqcut) \ +": nwt=" +sprinti("%3.0i", nwt) text = gsn_add_text(wks,plot(1),label,1850, 70, txres) plot(2) = gsn_csm_xy (wks,yyyy,ssnPha ,res) plot(3) = gsn_csm_xy (wks,yyyy,ssnPhau,res) ii = ind(.not.ismissing(ssnAmp)) ;******************************************** ; create attached plots ;******************************************** res1 = True res2 = True res2@gsnAttachPlotsXAxis = True amid = gsn_attach_plots(plot(0),(/plot(1),plot(2),plot(3)/),res1,res2) draw(plot) frame(wks)