WAFOFIG5 Joint distribution (pdf) of crest front velocity and wave height: Theoretical joint density of Vcf and 2*Ac (solid), kernel density estimate (dash) of Vcf and Hd data from Gullfaks C in the North Sea (dots) ( Increase NNp and Nh to get a smoother theoretical distribution. This may increase the computational time dramatically)
Estimate one-sided spectral density from data. | |
Extracts waveheights and steepnesses from data. | |
evaluates a PDF struct by interpolation | |
Binned Kernel Density Estimator. | |
Create or alter KDE OPTIONS structure. | |
Plot contents of pdf structures | |
Calculates quantile levels which encloses P% of PDF | |
Calculates quantile levels which encloses P% of data | |
Create or alter RIND OPTIONS structure. | |
Calculates a smoothing spline. | |
Joint density of amplitude and period/wave-length characteristics | |
Transforms process X and up to four derivatives | |
Prints a caption "made by WAFO" in current figure. | |
Clear variables and functions from memory. | |
String representation of date. | |
Find one string within another. | |
Hold current graph. | |
1-D interpolation (table lookup) | |
True if field is in structure array. | |
Linearly spaced vector. | |
X and Y arrays for 3-D plots. | |
Current date and time as date number. | |
Convert number to string. (Fast version) | |
Linear plot. | |
Display warning message; disable or enable warning messages. |
001 function wafofig5 002 % WAFOFIG5 Joint distribution (pdf) of crest front velocity and wave height: 003 % Theoretical joint density of Vcf and 2*Ac (solid), 004 % kernel density estimate (dash) of Vcf and Hd 005 % data from Gullfaks C in the North Sea (dots) 006 % 007 % ( Increase NNp and Nh to get a smoother theoretical distribution. 008 % This may increase the computational time dramatically) 009 010 011 % revised pab Feb2005 012 % -updated call to kdebin 013 014 global WAFOFIGNUM 015 016 if isempty(WAFOFIGNUM) 017 disp('You must start wafodemo in order to run this script') 018 clear global WAFOFIGNUM 019 return 020 end 021 022 global NVcf NHd Nxr Nrate 023 global fTcfAc NNp Nh Nnit Nspeed 024 global kdeVcfHd Nkernel Nhs NL2 025 026 if Nnit<0 027 opt = rindoptset('method',abs(Nnit),'speed',Nspeed); 028 else 029 opt = rindoptset('method',0,'nit',(Nnit),'speed',Nspeed); 030 end 031 % Only need to calculate Globals which is not empty 032 033 % Gullfaks C / North Sea data 034 % Extract Vcf and Hd 035 if isempty(NVcf) 036 [NVcf NHd] = dat2steep(Nxr,Nrate,0); 037 end 038 039 % Joint distribution of Tcf and Ac 040 if isempty(fTcfAc) 041 Sn=dat2spec(Nxr); 042 warning(['This takes several hours to finish ... (1 hour or more depending' ... 043 ' on input parameters and your computer)']) 044 fTcfAc = spec2thpdf(Sn,0,'TcfAc',[0 11 NNp],Nh,opt); 045 end 046 % Transform to the joint distribution of Vcf and 2*Ac 047 fVcf2Ac=th2vhpdf(fTcfAc); 048 049 050 % Kernel density estimate of Vcf and Hd 051 if isempty( kdeVcfHd) 052 kopt = kdeoptset('kernel',Nkernel,'hs',Nhs,'L2',NL2); 053 kdeVcfHd=kdebin([NVcf, NHd],kopt); 054 055 % calculate the levels which encloses fkde.pl percent of the data (v,h) 056 r = evalpdf(kdeVcfHd,NVcf, NHd,'linear'); 057 kdeVcfHd.cl = qlevels2(r,kdeVcfHd.pl); 058 end 059 060 061 plot( NVcf, NHd,'.'), hold on 062 pdfplot( kdeVcfHd,'r--') 063 pdfplot(fVcf2Ac,'k-') 064 hold off 065 wafostamp('Figure 5','(NR)') 066 067 return 068 069 function f2 = th2vhpdf(f,v) 070 % TH2VHPDF Calculates joint density of crest velocity and height 071 % in a stationary Gaussian transform process X(t) where 072 % Y(t) = g(X(t)) (Y zero-mean Gaussian with spectrum given in S). 073 % The transformation, g, can be estimated using lc2tr or dat2tr. 074 % 075 % CALL: fvh = th2vhpdf(fth,v); 076 % 077 % fth = (T Ac) pdf structure 078 % fvh = (V Ac)=(Ac/T Ac) pdf structure 079 % v = vector of velocities; note v >= 0, 080 % 081 % Example: 082 % fth = spec2thpdf(S,[],'TcfAc',[5 5 51],[],-2,4); 083 % fvh = th2vhpdf(fth); 084 % 085 % See also spec2thpdf 086 087 % Tested on : matlab 5.3 088 % History: 089 % revised by Per A. Brodtkorb 19.09.1999 090 091 %tic 092 093 f2=f; 094 f2.date=datestr(now); 095 f2.labx{1}='Velocity (m/s)'; 096 f2.title(find(f2.title=='T'))='V'; 097 k=findstr(f2.title,'min'); 098 if any(k) 099 f2.title=[f2.title(1:k(1)) 'ax' f2.title(k(1)+3:end)] ; 100 end 101 t=f.x{1}(:); 102 h=f.x{2}(:); 103 if nargin<2|isempty(v) 104 v=linspace(0,5,length(t))'; 105 else 106 v=v(:); 107 end 108 f2.x{1}=v; 109 110 for ix=2:length(h) 111 v1 =[0; h(ix)./t(2:end) ]; 112 %hold on,plot(v,f_th.f(ix,2:end).*t(2:end).^2/h(ix) ) 113 if 0 114 % the extrapolation done here does not work well if f.f si sparsely 115 % sampled => f2.f values corrupted 116 117 v1 =[0; h(ix)./t(2:end) ]; 118 ind= find(f.f(ix,:)>0); 119 f2.f(ix,2:end)=exp(min(smooth(v1(ind), log(f.f(ix,ind)'.*t(ind).^2/h(ix)) ,.99,v(2:end) ,1),0)); 120 else 121 % no extrapolation => more robust 122 v1 =[ h(ix)./t(2:end) ]; 123 124 f2.f(ix,2:end)=interp1(v1,[ f.f(ix,2:end)'.*t(2:end).^2/h(ix) ],v(2:end) ,'linear'); 125 end 126 end 127 f2.f(isnan(f2.f))=0; 128 k=find(f2.f>1000); 129 if any(k) 130 disp('Spurios spikes due to transformation') 131 %disp('set them to zero') 132 %f2.f(k)=0; 133 end 134 135 136 137 switch 2 138 case 0, %do nothing 139 case 2,% transform Amplitude 140 rate=2; 141 f2.f=f2.f/rate; 142 f2.x{2}=f2.x{2}*rate; 143 k=findstr(f2.title,'A'); 144 if any(k) 145 f2.title=[f2.title(1:k(1)-1) num2str(rate) f2.title(k(1):end)] ; 146 end 147 case 3 148 % this does not work correctly 149 g=f.tr; 150 utc=f.u; 151 152 Ac=f2.x{2}(:); 153 Nx=length(Ac); 154 Nv=length(f2.x{1}(:)); 155 der=ones(Nx,1); % dh/dh=1 156 hg=tranproc([utc+Ac der],g); 157 der1=hg(:,2); 158 Acg=hg(:,1); % Gaussian level 159 hg=tranproc([-Acg der],fliplr(g)); 160 der2=hg(:,2); 161 der=1+abs(der1.*der2); 162 % transform amplitude to wave height H=Ac-G(-g(Ac+utc)); 163 f2.x{2}=Ac-tranproc(-Acg,fliplr(g)); 164 f2.f=f2.f.*der(:,ones(1,size(f2.f,2))); 165 x1=linspace(f2.x{1}(1),f2.x{1}(Nv),Nv)'; 166 x2=linspace(f2.x{2}(1),f2.x{2}(Nx),Nx)'; 167 [X1 X2]=meshgrid(x1,x2); 168 f2.f = evalpdf(f2,X1,X2,'linear'); 169 f2.x{2}=x2; 170 f2.x{1}=x1; 171 172 k=findstr(f2.title,'Ac'); 173 if any(k) 174 f2.title=[f2.title(1:k(1)-1) 'Ac-G(-g(Ac+utc))' f2.title(k(1)+2:end)] ; 175 end 176 end 177 178 if isfield(f,'pl') 179 f2.cl=qlevels(f2.f,f.pl); 180 else 181 f2.pl=[10:20:90 95 99 99.9]; 182 f2.cl=qlevels(f2.f,f2.pl); 183 end 184 %toc
Comments or corrections to the WAFO group