0001 function [f] = spec2Acdf(spec,utc,def,paramt,h,nit,speed,bound)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053
0054
0055
0056
0057 startTime = clock;
0058 if nargin<3|isempty(def)
0059 def='tc';
0060 end
0061 if strcmp('l',lower(def(1)))
0062 spec=spec2spec(spec,'k1d');
0063 elseif strcmp('t',lower(def(1)))
0064 spec=spec2spec(spec,'freq');
0065 else
0066 error('Unknown def')
0067 end
0068 switch lower(def)
0069 case {'tc','lc'}, defnr = 1;
0070 case {'tt','lt'}, defnr =-1;
0071 otherwise, ,error('Unknown def')
0072 end
0073
0074 [S, xl4, L0, L2, L4, L1]=wnormspec(spec);
0075
0076 A = sqrt(L0/L2);
0077 SCIS=0;
0078 if nargin<6|isempty(nit)
0079 nit=2;
0080 elseif nit<0
0081 SCIS=min(abs(nit),2);
0082 nit=1;
0083 end
0084
0085 if isfield(spec,'tr')
0086 g=spec.tr;
0087 else
0088 g=[];
0089 end
0090 if isempty(g)
0091 g=[sqrt(L0)*(-5:0.02:5)', (-5:0.02:5)'];
0092 end
0093
0094
0095
0096 if nargin<2|isempty(utc)
0097 utc_d = gaus2dat([0, 0],g);
0098 utc = utc_d(1,2);
0099 end
0100
0101
0102 uu = dat2gaus([0., utc],g);
0103 u = uu(2);
0104 disp(['The level u for Gaussian process = ', num2str(u)])
0105
0106
0107
0108 if nargin<7|isempty(speed)
0109 speed=4;
0110 end
0111
0112 if nargin<8|isempty(bound)
0113 bound=0;
0114 end
0115
0116 if SCIS>0
0117 bound=0;
0118 end
0119 ttn = 1.2*ceil(2*pi*sqrt(L0/L2)*exp(u^2/2)*(0.5+erf(-sign(defnr)*u/sqrt(2))/2));
0120 if nargin<4|isempty(paramt)
0121 Ntime = 51;
0122 t0 = ttn;
0123 tn = 2*ceil(2*pi*sqrt(L0/L2)*exp(u^2/2)*(0.5+erf(-sign(defnr)*u/sqrt(2))/2));
0124 paramt=[t0,tn,Ntime];
0125 else
0126 t0 = paramt(1);
0127 tn = paramt(2);
0128 Ntime = paramt(3);
0129 end
0130
0131 t0=ttn;
0132 t0=ttn/A;
0133 tt = tn/A;
0134 t = levels([0, tt, Ntime]);
0135 N0 = 2+ceil(ttn/tn*(Ntime-1));
0136 NN0=N0;
0137 if Ntime>101
0138 disp('nr. of wavelengths > 101., slow')
0139 end
0140
0141
0142 nr = 2;
0143
0144 if nargin<5|isempty(h)
0145 hg=[0:0.25:5.];
0146 h=tranproc(u+sign(defnr)*hg',fliplr(g))-utc;
0147 end
0148
0149
0150 h=reshape(h,length(h),1);
0151 Nx=length(h);
0152 if Nx>101
0153 error('nr. of amplitudes limited to 101.')
0154 end
0155 hg=tranproc([utc+sign(defnr)*h],g);
0156
0157 if exist('h.in'), delete h.in, end
0158 fid=fopen('h.in','wt');
0159 fprintf(fid,'%12.10f\n',hg);
0160 fclose(fid);
0161
0162
0163 dt = t(2)-t(1);
0164
0165 R = spec2cov2(S,nr,Ntime-1,dt);
0166
0167
0168
0169
0170
0171
0172
0173
0174
0175
0176
0177 for k=0:nr
0178 filename=['Cd' int2str(k) '.in'];
0179 if exist(filename)
0180 delete(filename)
0181 end
0182 fid=fopen(filename,'wt');
0183 fprintf(fid,'%12.10E \n',R(:,k+1));
0184 fclose(fid);
0185 end
0186
0187
0188 if exist('reflev.in'), delete reflev.in, end
0189 disp('writing data')
0190 fid=fopen('reflev.in','wt');
0191 fprintf(fid,'%12.10E \n',u);
0192 defnr;
0193
0194 fprintf(fid,'%2.0f \n',defnr);
0195 fprintf(fid,'%2.0f \n',Ntime);
0196 fprintf(fid,'%2.0f \n',N0);
0197 fprintf(fid,'%2.0f \n',nit);
0198 fprintf(fid,'%2.0f \n',speed);
0199 fprintf(fid,'%2.0f \n',SCIS);
0200 fprintf(fid,'%2.0f \n',10^9*rand);
0201 fprintf(fid,'%2.0f \n',Nx);
0202 fprintf(fid,'%12.10E \n',dt);
0203 fclose(fid);
0204 disp(' Starting Fortran executable.')
0205
0206 if bound<0.5
0207 dos([ wafoexepath, 'sp2Acdf70.exe']);
0208 else
0209 dos([ wafoexepath, 'sp2Acdf51.exe']);
0210 end
0211
0212
0213 switch defnr
0214 case 1, Htxt = ['Distribution of crest height'];
0215 case -1, Htxt = ['Distribution of trough height'];
0216 end
0217
0218 f = createpdf;
0219 f.title=Htxt;
0220 f.nit=nit;
0221 f.speed=speed;
0222 f.SCIS=SCIS;
0223 f.u=utc;
0224 f.labx1='amplitude [m]';
0225 ftmp=loaddata('dens.out')/A;
0226 f.x{1}=h;
0227
0228
0229
0230 if bound<0.5
0231 ftmp=reshape(ftmp,Nx,length(t));
0232 if (length(t)>2)
0233 ftmp(:,2)=0.5*(ftmp(:,3)+ftmp(:,1));
0234 end
0235 else
0236 ftmp_up=reshape(ftmp(:,1),Nx,length(t));
0237 ftmp_lo=reshape(ftmp(:,2),Nx,length(t));
0238 if (length(t)>2)
0239 ftmp_up(:,2)=0.5*(ftmp_up(:,3)+ftmp_up(:,1));
0240 ftmp_lo(:,2)=0.5*(ftmp_lo(:,3)+ftmp_lo(:,1));
0241 end
0242 end
0243
0244
0245
0246
0247
0248
0249
0250
0251 Ntime = 61;
0252 t0 = 0.;
0253 tn = ttn;
0254 paramt1=[t0, tn, Ntime];
0255 tt = tn/A;
0256 t = levels([0, tt, Ntime]);
0257 N0 = 1+ceil(t0/tn*(Ntime-1));
0258 dt=t(2)-t(1);
0259 R = spec2cov2(S,nr,Ntime-1,dt);
0260 for k=0:nr
0261 filename=['Cd', int2str(k), '.in'];
0262 if exist(filename)
0263 delete(filename)
0264 end
0265 fid = fopen(filename,'wt');
0266 fprintf(fid,'%12.10E \n',R(:,k+1));
0267 fclose(fid);
0268 end
0269 if exist('reflev.in'), delete reflev.in, end
0270 disp('writing data')
0271 fid=fopen('reflev.in','wt');
0272 fprintf(fid,'%12.10E \n',u);
0273 defnr;
0274
0275 fprintf(fid,'%2.0f \n',defnr);
0276 fprintf(fid,'%2.0f \n',Ntime);
0277 fprintf(fid,'%2.0f \n',N0);
0278 fprintf(fid,'%2.0f \n',nit);
0279 fprintf(fid,'%2.0f \n',speed);
0280 fprintf(fid,'%2.0f \n',SCIS);
0281 fprintf(fid,'%2.0f \n',10^9*rand);
0282 fprintf(fid,'%2.0f \n',Nx);
0283 fprintf(fid,'%12.10E \n',dt);
0284 fclose(fid);
0285 disp(' Starting Fortran executable.')
0286 if bound<0.5
0287 dos([ wafoexepath, 'sp2Acdf70.exe']);
0288 else
0289 dos([ wafoexepath, 'sp2Acdf51.exe']);
0290 end
0291
0292
0293 ftmp1=loaddata('dens.out')/A;
0294
0295 if bound<0.5
0296 ftmp1=reshape(ftmp1,Nx,length(t));
0297 if (length(t)>2)
0298 ftmp1(:,2)=0.5*(ftmp1(:,3)+ftmp1(:,1));
0299 end
0300 else
0301 ftmp_up1=reshape(ftmp1(:,1),Nx,length(t));
0302 ftmp_lo1=reshape(ftmp1(:,2),Nx,length(t));
0303 if (length(t)>2)
0304 ftmp_up1(:,2)=0.5*(ftmp_up1(:,3)+ftmp_up1(:,1));
0305 ftmp_lo1(:,2)=0.5*(ftmp_lo1(:,3)+ftmp_lo1(:,1));
0306 end
0307 end
0308
0309
0310 X=[levels(paramt1)];
0311 Nt=paramt(3);
0312 if (NN0<Nt+1)
0313 t=levels([0, paramt(2), paramt(3)]);
0314 X=[t(NN0:Nt), X];
0315 N=length(X);
0316 [X, indX]=sort(X);
0317 dX=X;
0318 dx(1)=0.;
0319 dx(N)=X(N)-X(N-1);
0320 dx(2:N-1)=0.5*(X(3:N)-X(1:N-2));
0321 if bound<0.5
0322 Y=[ftmp(:,NN0:Nt), ftmp1];
0323 Y=Y(:,indX);
0324 f.f=min(Y*dx',1);
0325
0326
0327 else
0328 Y_up=[ftmp_up(:,NN0:Nt), ftmp_up1];
0329 Y_up=Y_up(:,indX);
0330 Y_lo=[ftmp_lo(:,NN0:Nt), ftmp_lo1];
0331 Y_lo=Y_lo(:,indX);
0332 f.f=min([Y_up*dx', Y_lo*dx'],1);
0333
0334
0335 end
0336 else
0337 N=length(X);
0338 dX=X;
0339 dx(1)=0.;
0340 dx(N)=X(N)-X(N-1);
0341 dx(2:N-1)=0.5*(X(3:N)-X(1:N-2));
0342 if bound<0.5
0343 Y=[ftmp1];
0344 f.f=min(Y*dx',1);
0345
0346
0347 else
0348 Y_up=[ftmp_up1];
0349 Y_lo=[ftmp_lo1];
0350 f.f=min([Y_up*dx', Y_lo*dx'],1);
0351 end
0352 end
0353 f.elapsedTime = etime(clock,startTime);
0354
0355 pdfplot(f)
0356 hold on
0357 if abs(u)<0.1
0358 xu=tranproc(u+f.x{1},g);
0359 plot(f.x{1},1-exp(-xu.*xu/2),'b.')
0360 plot(f.x{1},1-exp(-xu.*xu/2))
0361 end
0362