001 function recinit
002
003
004
005
006
007
008
009
010
011
012
013
014
015
016
017
018
019
020
021
022
023
024
025
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;
037 end
038
039 if isempty(recmenulabels)
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)
053 xn=load(recfilename);
054
055 end
056 dT=xn(2,1)-xn(1,1);
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
068
069 global inds zcrit dcrit ddcrit
070 global Nsim L csm1 csm2 param
071 global xr g g2 test0 tobs mu1o mu1oStd
072 global Sr m0 m2 Tm02 Hm0 Vrms Hrms
073 global V H rate
074 global phat res noverlap
075 global sphat CSMA CSMB
076 global ft2 fkde
077 global kernel hs L2
078
079
080
081
082
083
084
085 if isempty(zcrit), zcrit = 0.02; end
086
087 if isempty(dcrit), dcrit = 5*dT; end
088 if isempty(ddcrit), ddcrit = gravity(61)*dT^2; end
089
090
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
101
102
103
104
105 if isempty(rate), rate = 8; end
106
107
108
109
110 if isempty(res), res = 0.2; end
111 if isempty(noverlap); noverlap = 0; end
112
113
114
115 if isempty(CSMA), CSMA = 0.95; end
116 if isempty(CSMB), CSMB = 0.95; end
117
118
119
120 if isempty(kernel), kernel = 'epan'; end
121 if isempty(L2), L2 = [.5 .5]; end
122
123
124
125
126
127
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
142 Sr=dat2spec(xr,L);
143 m0=trapz(Sr.w,Sr.S);
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
154
155
156
157 if RECFIGNUM>6,
158 if isempty(V) | isempty(H)
159
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 ]);
173 end
174
175
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,
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);
189
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);
198
199 end
200 fkde.labx{1}='v';
201 fkde.labx{2}='h';
202 end
203 end