Home > wafo > papers > rec > recinit.m

recinit

PURPOSE ^

setup all global variables of the RECDEMO

SYNOPSIS ^

recinit

DESCRIPTION ^

 RECINIT  setup all global variables of the RECDEMO

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function recinit
002 %RECINIT  setup all global variables of the RECDEMO
003 
004 % Since the calculations are very dependent on eachother,
005 % all the time consuming calculations are located here.
006 % Plotting and some minor calculations are left to the recfigXX.m
007 % functions 
008 %
009 % GLOBAL RECFIGNUM <=2 : globals defined and datasets are loaded
010 %                  >=3 : reconstruction of xn, extraction of V and H and
011 %                        distribution fitting to V and H
012 %
013 % Recinit only recalculates the values of empty 
014 % variables if RECFIGNUM>=3 ==> saves alot of computation time
015 
016 % revised pab feb2005
017 % updated call to kdebin
018 % revised pab feb2004
019 % revised pab 18.04.2001
020 % - updated call to reconstruct due to changes in all transformation
021 %    estimation functions
022 % - updated call to dat2steep
023 % By pab 28.01.2000
024 
025 % Figure parameter:
026 %~~~~~~~~~~~~~~~~~~
027 global RECFIGNUM  
028 if isempty(RECFIGNUM)
029   disp('You must start recdemo in order to run this script')
030   clear global RECFIGNUM
031   return
032 end
033 
034 global pwdstr  recmenulabels recfilename xn map wind 
035 if isempty(pwdstr)
036   pwdstr=pwd; % save path to where the demo was started
037 end
038 
039 if isempty(recmenulabels) % automatic generation of menu labels
040   cd(fullfile(waforoot,'papers','rec'));  
041   recmenulabels=cell(12,1);
042   for ix=1:13
043     recmenulabels{ix} = geth1line(['recfig' num2str(ix)],1);   
044   end
045   cd(pwdstr)
046 end
047 
048 if isempty(recfilename),
049   recfilename='gfaks89.dat';
050 end
051 
052 if isempty(xn) % load if not already loaded 
053   xn=load(recfilename);
054   %xn=xn(1:10000,:); % used for debugging 
055 end
056 dT=xn(2,1)-xn(1,1); % sampling
057 
058 if isempty(map)
059   disp('Loading map over the North Sea...')
060   map = load('northsea.dat');
061 end
062 if isempty(wind)
063   disp('Loading wind measurements from Statfjord A...')
064   wind = load('sfa89.dat');
065 end
066 
067 % Define Globals used
068 %~~~~~~~~~~~~~~~~~~~~~~~~~~~
069 global inds zcrit dcrit ddcrit        % findoutliers output, input
070 global Nsim L csm1 csm2 param         % reconstruct input
071 global xr g g2 test0 tobs mu1o mu1oStd % reconstruct output                
072 global Sr m0 m2 Tm02 Hm0 Vrms Hrms    % Spectral density, moments and Sea state parameters
073 global V H rate                       % dat2steep output, input
074 global phat res noverlap              % dist2dfit output, input
075 global sphat CSMA CSMB                % dist2dsmfun2 output, input
076 global ft2 fkde                       % dist2dpdf2 and kdebin outputs
077 global kernel hs L2                   % kdebin input: kernel, smoothing parameter and transformation, respectively
078 
079 
080 % RECONSTRUCTION 
081 %~~~~~~~~~~~~~~~~
082 
083 % Set findoutliers default values:
084 %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
085 if isempty(zcrit),  zcrit  = 0.02;               end
086 % set dcrit corresponding to 5m/s and ddcrit corresponding to g (maximum acceleration of Stokes wave is g/2) 
087 if isempty(dcrit),  dcrit  = 5*dT;               end % dcrit=1.23 used in  article
088 if isempty(ddcrit), ddcrit = gravity(61)*dT^2;   end % ddcrit=1.5*1.23 used in  article
089 
090 % Set reconstruct default values:
091 %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
092 if isempty(Nsim),   Nsim   = 6;                  end 
093 if isempty(csm1),   csm1   = 0.95;               end
094 if isempty(csm2),   csm2   = 0.05;               end
095 if isempty(param),  param  = [-5 5 501];        end
096 
097 
098 %  
099 
100 % DISTRIBUTION FITTING
101 %~~~~~~~~~~~~~~~~~~~~~
102 
103 % set dat2steep default value
104 %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
105 if isempty(rate),   rate   = 8;                  end % interpolation rate of data set
106                                                      %  before extraction of parameters
107 %
108 %Set dist2dfit default values:
109 %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
110 if isempty(res),        res    = 0.2;  end
111 if isempty(noverlap); noverlap = 0;    end
112 
113 %Set dist2dsmfun2 default values:
114 %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
115 if isempty(CSMA),         CSMA = 0.95; end
116 if isempty(CSMB),         CSMB = 0.95; end
117                              
118 % Set kdebin default values: 
119 %~~~~~~~~~~~~~~~~~~~~~~~~~~~
120 if isempty(kernel),     kernel = 'epan';  end
121 if isempty(L2),             L2 = [.5 .5]; end
122 % if hs and L2 is empty then kdebin uses the 
123 % default values calculated  from the data (v,h)                             
124 
125 
126 
127 % RECONSTRUCTION 
128 %~~~~~~~~~~~~~~~~
129 
130 if RECFIGNUM>2,
131   if isempty(inds)
132     disp('Finding spurious points')
133     inds = findoutliers(xn,zcrit,dcrit,ddcrit);
134   end
135   if isempty(xr)
136     disp('Reconstruction ...')
137     disp('This may take a while (20 min -> 1 hour depending on machine, inputs etc...)')
138     [xr,g,g2,test0,tobs,mu1o,mu1oStd]=reconstruct(xn,inds,Nsim,L,[],...
139     troptset('csm',csm1,'gsm',csm2,'param',param));
140   
141     % Estimate spectral density and moments
142     Sr=dat2spec(xr,L);
143     m0=trapz(Sr.w,Sr.S); % or alternatively spec2mom
144     m2=trapz(Sr.w,Sr.S);
145     Hm0=4*sqrt(m0);
146     Tm02=2*pi*sqrt(m0/m2);
147     Vrms=2*Hm0/Tm02;
148     Hrms=Hm0/sqrt(2);
149   end
150 end
151 
152 
153 % DISTRIBUTION FITTING:
154 %~~~~~~~~~~~~~~~~~~~~~~
155 
156 
157 if RECFIGNUM>6,                 
158   if isempty(V) | isempty(H)
159     % Not extracting the missing points   
160     Ns=min(find(isnan(xn(:,2)))) 
161     if ~isempty(Ns), [V , H]=dat2steep(xr(1:Ns-1,:),rate,0);end
162     Ne=max(find(isnan(xn(:,2)))); 
163     if isempty(Ne),Ne=0;end
164     [V2 , H2]=dat2steep(xr(Ne+1:end,:),rate,0);
165     V=[V(:);V2(:)];H=[H(:);H2(:)];
166   end
167   if isempty(phat)
168     disp('Fitting a 2D distribution to the data')
169     phat=dist2dfit(V/Vrms,H/Hrms,{'weibull','rayleigh'},[res,noverlap]);
170   end
171   if isempty(sphat)
172     sphat=dist2dsmfun2(phat,linspace(0,3.5),[CSMA,CSMB ],[1 1 ]); % smooth the estimated parameter values
173   end
174 
175   % Calculate the theoretical density of fitted model
176   if isempty(ft2)                    
177     ft2=dist2dpdf2(linspace(0,4),linspace(0,4),sphat);
178     ft2.labx{1}='v';
179     ft2.labx{2}='h';
180   end
181 end
182 
183 
184 if RECFIGNUM>8, % kernel density estimation
185   if isempty(fkde)
186     disp('Estimating kernel density estimate')
187     kopt = kdeoptset('kernel',kernel,'L2',L2,'inc',128);
188     fkde = kdebin([V/Vrms,H/Hrms],kopt); % kdebin is much
189                                                         % faster than kde
190     k=find(fkde.f>30);
191     if any(k)
192       disp('Removing the spurios spikes')
193       fkde.f(k)=0;
194     end
195     if 1,
196       r = evalpdf(fkde,V/Vrms,H/Hrms,'linear');
197       fkde.cl = qlevels2(r,fkde.pl); % calculate the levels which encloses fkde.pl
198       % percent of the data (v,h)
199     end
200     fkde.labx{1}='v';
201     fkde.labx{2}='h';
202   end
203 end

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