TH2VHPDF Transform joint T-H density to V-H density CALL: fvh = th2vhpdf(fth,v); fth = (T Ac) pdf structure fvh = (V Ac)=(Ac/T Ac) pdf structure v = vector of velocities; note v >= 0, Example: opt = rindoptset('speed',5,'nit',1,'method',0); S = jonswap; fth = spec2thpdf(S,4,'TcfAc',[0 5 51],[],opt); fvh = th2vhpdf(fth); See also spec2thpdf
evaluates a PDF struct by interpolation | |
Calculates quantile levels which encloses P% of PDF | |
Calculates a smoothing spline. | |
Transforms process X and up to four derivatives | |
String representation of date. | |
Find one string within another. | |
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) |
001 function f2 = th2vhpdf(f,v) 002 %TH2VHPDF Transform joint T-H density to V-H density 003 % 004 % CALL: fvh = th2vhpdf(fth,v); 005 % 006 % fth = (T Ac) pdf structure 007 % fvh = (V Ac)=(Ac/T Ac) pdf structure 008 % v = vector of velocities; note v >= 0, 009 % 010 % Example: 011 % opt = rindoptset('speed',5,'nit',1,'method',0); 012 % S = jonswap; 013 % fth = spec2thpdf(S,4,'TcfAc',[0 5 51],[],opt); 014 % fvh = th2vhpdf(fth); 015 % 016 % See also spec2thpdf 017 018 % Tested on : matlab 5.3 019 % History: 020 % revised pab 3Dec2003 021 % revised by Per A. Brodtkorb 19.09.1999 022 023 %tic 024 025 f2=f; 026 f2.date=datestr(now); 027 f2.labx{1}='Velocity (m/s)'; 028 f2.title(find(f2.title=='T'))='V'; 029 k=findstr(f2.title,'min'); 030 if any(k) 031 f2.title=[f2.title(1:k(1)) 'ax' f2.title(k(1)+3:end)] ; 032 end 033 t=f.x{1}(:); 034 h=f.x{2}(:); 035 if nargin<2|isempty(v) 036 v=linspace(0,4,length(t))'; 037 else 038 v=v(:); 039 end 040 f2.x{1}=v; 041 f2.f = zeros(length(f2.x{2}),length(v)); 042 k0 = find(v); 043 k = find(t); 044 k3 = find(h(:)).'; 045 method = 'linear' 046 for ix=k3 047 if 0 048 % the extrapolation done here does not work well if f.f si sparsely 049 % sampled => f2.f values corrupted 050 051 v1 = h(ix)./t(k) ; 052 ind = find(f.f(ix,k)>0); 053 f2.f(ix,k0)=exp(min(smooth(v1(ind), log(f.f(ix,k(ind)).'.*t(k(ind)).^2/h(ix)) ,.99,v(k0) ,1),0)); 054 else 055 % no extrapolation => more robust 056 v1 =[ h(ix)./t(k) ]; 057 f2.f(ix,k0)=interp1(v1,f.f(ix,k).'.*t(k)./v1 ,v(k0) ,method); 058 end 059 end 060 f2.f(isnan(f2.f))=0; 061 k=find(f2.f>1000); 062 if any(k) 063 disp('Spurios spikes due to transformation') 064 %disp('set them to zero') 065 %f2.f(k)=0; 066 end 067 068 069 070 switch 0 071 case 0, %do nothing 072 case 2,% transform Amplitude 073 rate=1.9 074 f2.f=f2.f/rate; 075 f2.x{2}=f2.x{2}*rate; 076 k=findstr(f2.title,'A'); 077 if any(k) 078 f2.title=[f2.title(1:k(1)-1) num2str(rate) f2.title(k(1):end)] ; 079 end 080 case 3 081 g=f.tr; 082 utc=f.u; 083 084 Ac=f2.x{2}(:); 085 Nx=length(Ac); 086 Nv=length(f2.x{1}(:)); 087 der=ones(Nx,1); % dh/dh=1 088 hg=tranproc([utc+Ac der],g); 089 der1=hg(:,2); 090 Acg=hg(:,1); % Gaussian level 091 hg=tranproc([-Acg der],fliplr(g)); 092 der2=hg(:,2); 093 der=1+abs(der1.*der2); 094 % transform amplitude to wave height H=Ac-G(-g(Ac+utc)); 095 f2.x{2}=Ac-tranproc(-Acg,fliplr(g)); 096 f2.f=f2.f.*der(:,ones(1,size(f2.f,2))); 097 x1=linspace(f2.x{1}(1),f2.x{1}(Nv),Nv)'; 098 x2=linspace(f2.x{2}(1),f2.x{2}(Nx),Nx)'; 099 [X1 X2]=meshgrid(x1,x2); 100 f2.f = evalpdf(f2,X1,X2,'linear'); 101 f2.x{2}=x2; 102 f2.x{1}=x1; 103 104 k=findstr(f2.title,'Ac'); 105 if any(k) 106 f2.title=[f2.title(1:k(1)-1) 'Ac-G(-g(Ac+utc))' f2.title(k(1)+2:end)] ; 107 end 108 end 109 110 if isfield(f,'pl') 111 f2.cl=qlevels(f2.f,f.pl); 112 else 113 f2.pl=[10:20:90 95 99 99.9]; 114 f2.cl=qlevels(f2.f,f2.pl); 115 end 116 %toc
Comments or corrections to the WAFO group