Home > wafo > trgauss > th2vhpdf.m

# th2vhpdf

## PURPOSE

Transform joint T-H density to V-H density

## SYNOPSIS

f2 = th2vhpdf(f,v)

## DESCRIPTION

``` 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);

## CROSS-REFERENCE INFORMATION

This function calls:
 evalpdf evaluates a PDF struct by interpolation qlevels Calculates quantile levels which encloses P% of PDF smooth Calculates a smoothing spline. tranproc Transforms process X and up to four derivatives datestr String representation of date. findstr Find one string within another. interp1 1-D interpolation (table lookup) isfield True if field is in structure array. linspace Linearly spaced vector. meshgrid X and Y arrays for 3-D plots. now Current date and time as date number. num2str Convert number to string. (Fast version)
This function is called by:

## SOURCE CODE

```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 %
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```

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