Dear Ferret Users,
I
need further suggestions to solve my problem in overlaying of the
error_bars or plot_swath into the power spectrum of the time series observation and the
simulation.
Since I need the horizontal axes in period and vertical axes in log scale, I couldn't find the right script.
I always got the results of the error_bars or plot_swath in the meter (for vertical axis) or frequency (for the horizontal axis).
Hopefully there's a way to solve these case.
Following is the script which the grey font shows that I'm not sure with those script:
!script for power spectrum and error_bars
!-----
cancel mode logo
use
"/media/sep052012/verification.cdf"
let FFT_simulation=simulation
let FFT_simulationfft=FFTA(FFT_simulation[l=1897:2640])
set variable/title="Amplitude Spectrum" FFT_simulationfft
let FFT_verification=verification
let FFT_verificationfft=FFTA(FFT_verification[l=1:744])
set variable/title="Amplitude Spectrum" FFT_verificationfft
!==simulation
LET FFT_nf = `FFT_simulationfft,return=lend`
LET FFT_nyquist = 0.5
LET FFT_freq1 = FFT_nyquist/ FFT_nf
DEFINE AXIS/T=`FFT_freq1`:`FFT_nyquist`:`FFT_freq1` faxis
DEFINE GRID/T=faxis gfftfreq
LET a = T[g=gfftfreq]
Let per = 1./a
!==verification
LET FFT_nf2 = `FFT_verificationfft,return=lend`
LET FFT_nyquist2 = 0.5
LET FFT_freq2 = FFT_nyquist2/ FFT_nf2
DEFINE AXIS/T=`FFT_freq2`:`FFT_nyquist2`:`FFT_freq2` faxis2
DEFINE GRID/T=faxis2 gfftfreq2
LET b = T[g=gfftfreq2]
Let per2 = 1./b
!==plot power spectrum in log-log scale
set win 1
set win/size=1./aspect=1.:axis
PLOT/vlog/VS/dash/line=1/vlimits=0.00001:10/HLIMITS=0:40:2/thick=3/nolab per[l=1:`FFT_nf`], FFT_simulationfft
PPL YLAB Amplitude (m) log scale
ppl plot
PLOT/vlog/OVERLAY/VS/LINE=1/thick=3/vlimits=0.00001:10/HLIMITS=0:40:2/nolab per2[l=1:`FFT_nf2`], FFT_verificationfft
ppl color 6 100 100 100
label/nou -2.5 10.0 .03 "Tidal amplitude spectrum"
label/nou -2 9 .03 M_2
label/nou -2.5 8.8 .03 S_2
label/nou 0.5 8.9 .03 K_1
label/nou 1. 8.8 .03 O_1
label/nou -0.2 6.2 .03 PERIOD (HOURS)
!plot overlay error_bars and/or plot_swath
set win 2
!let/quiet errup = FFT_verificationfft + FFT_simulationfft
!let/quiet errdn =
FFT_verificationfft - FFT_simulationfft
let/quiet err = ((FFT_verificationfft[l=1:`FFT_nf2`] - FFT_simulationfft[l=1:`FFT_nf`])/2)
plot/vlog/VS/line=1/vlimits=0.00001:10/HLIMITS=0:40:2/thick=3/nolab per2[l=1:`FFT_nf2`], FFT_verificationfft
ppl xlab period (hours)
PPL YLAB Amplitude (m) log scale
ppl plot
!go plot_swath poly/ov/pal=red/nolab FFT_verificationfft-err FFT_verificationfft+err
!plot/vlog/vs/over/nolab/color=red/vlimits=0.00001:10/HLIMITS=0:40:2/thickness=2 per3[l=1:`FFT_nf3`], err
set win 3
!plot/vlog/VS/line=1/vlimits=0.00001:10/HLIMITS=0:40:2/thick=3/nolab FFT_verificationfft, err
!ppl xlab period (hours)
!PPL YLAB Amplitude (m^2) log scale
!ppl plot
!go error_bars poly/color=red/nolab FFT_verificationfft err
LET FFT_nf3 = `err,return=lend`
LET FFT_nyquist3 = 0.5
LET FFT_freq3 = FFT_nyquist3/ FFT_nf3
DEFINE AXIS/T=`FFT_freq3`:`FFT_nyquist3`:`FFT_freq3` faxis3
DEFINE
GRID/T=faxis3 gfftfreq3
LET c = T[g=gfftfreq3]
Let per3 = 1./c
!plot/vlog/VS/line=1/vlimits=0.000001:10/HLIMITS=0:40:2/thick=3/nolab per2[l=1:`FFT_nf2`], FFT_verificationfft
!go error_bars poly/ov/color=red/nolab FFT_verificationfft FFT_simulationfft
!ppl xlab period (hours)
!PPL YLAB Amplitude (m) log scale
!ppl plot
!plot/vlog/ov/VS/line=1/dash/vlimits=0.000001:10/HLIMITS=0:40:2/thick=3/nolab per[l=1:`FFT_nf`], FFT_simulationfft
Thank you in advance for any further suggestions to solve such problem.
Best wishes,
dessy