Home > wafo > papers > wafodemo > wafofig5.m

wafofig5

PURPOSE ^

Joint distribution (pdf) of crest front velocity and wave height:

SYNOPSIS ^

wafofig5

DESCRIPTION ^

  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)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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

Mathematical Statistics
Centre for Mathematical Sciences
Lund University with Lund Institute of Technology

Comments or corrections to the WAFO group


Generated on Thu 06-Oct-2005 02:21:16 for WAFO by m2html © 2003