Home > wafo > wsim > lc2sdat.m

lc2sdat

PURPOSE ^

Simulates process with given irregularity factor and crossing spectrum

SYNOPSIS ^

process=lc2sdat(lc,N,alpha)

DESCRIPTION ^

 LC2SDAT Simulates process with given irregularity factor and crossing spectrum 
 
   CALL: L = lc2sdat(lc,N,alpha);
 
         L     = two column matrix with times and values of the simulated
                 process.
         lc    = two column matrix with levels and crossing intensities,
                 i.e., level-crossing spectrum.. 
         N     = number of sample points.
         alpha = irregularity factor, 0<alpha<1, small  alpha  gives 
                 irregular process.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function process=lc2sdat(lc,N,alpha)
002 %LC2SDAT Simulates process with given irregularity factor and crossing spectrum 
003 %
004 %  CALL: L = lc2sdat(lc,N,alpha);
005 %
006 %        L     = two column matrix with times and values of the simulated
007 %                process.
008 %        lc    = two column matrix with levels and crossing intensities,
009 %                i.e., level-crossing spectrum.. 
010 %        N     = number of sample points.
011 %        alpha = irregularity factor, 0<alpha<1, small  alpha  gives 
012 %                irregular process.
013 %
014 
015 % Example
016 %  n = 10000;  
017 %  S = jonswap(7);
018 %  alpha = spec2char(S,'alpha');
019 %  xs  = spec2sdat(S,n);  
020 %  lc  = dat2lc(xs);
021 %  xs2 = lc2sdat(lc,n,alpha);
022 %  Se  = dat2spec(xs2);  
023 %  wspecplot(S),hold on
024 %  wspecplot(Se,'r')  
025   
026 %History
027 % revised pab Feb2004
028 % -added example  
029 % -changed order of input to conform with the other XXX2sdat routines  
030 % Revised by GL, July 10, 2000:
031 % Changed call, from  cross2tr  to  lc2tr
032 % Changed "Check the result" by removing reference to getrfc
033 
034 % TODO % add a good example
035   
036 error(nargchk(1,3,nargin))
037 
038 f = linspace(0,0.49999,1000);
039 rho_st = 2*sin(f*pi).^2-1;
040 tmp    = alpha*asin(sqrt((1+rho_st)/2));
041 tmp    = sin(tmp).^2;
042 a2     = (tmp-rho_st)./(1-tmp);
043 y      = min([a2+rho_st;1-a2]);
044 [maximum,maxidx]=max(y);
045 
046 rho_st = rho_st(maxidx);
047 a2     = a2(maxidx);
048 a1     = 2*rho_st+a2-1;
049 r0     = 1;
050 r1     = -a1/(1+a2);
051 r2     = (a1^2-a2-a2^2)/(1+a2);
052 sigma2 = r0+a1*r1+a2*r2;
053 
054 e      = randn(N,1)*sqrt(sigma2);
055 e(1:2) = [0;0];
056 L0     = randn(1,1);
057 L0     = [L0;r1*L0+sqrt(1-r2^2)*randn(1,1)];
058 %Simulate the process, starting in L0
059 L      = filter(1,[1 a1 a2],e,filter([1 a1 a2],1,L0));
060 
061 epsilon = 1.01;
062 min_L   = min(L);
063 max_L   = max(L);
064 maxi    = max(abs([min_L max_L]))*epsilon;
065 mini    = -maxi;
066 
067 u = linspace(mini,maxi,101)';
068 G = (1+erf(u/sqrt(2)))/2;
069 G = G.*(1-G);
070 
071 x = linspace(0,r1,100)';
072 factor1  = 1./sqrt(1-x.^2);
073 factor2  = 1./(1+x);
074 integral = zeros(size(u));
075 for i=1:length(integral)
076   y = factor1.*exp(-u(i)*u(i).*factor2);
077   integral(i) = trapz(x,y);
078 end
079 G = G-integral/(2*pi);
080 G = G/max(G);
081 
082 Z = ((u>=0)*2-1).*sqrt(-2*log(G));
083 
084 sumcr   = trapz(lc(:,1),lc(:,2));
085 lc(:,2) = lc(:,2)/sumcr;
086 mcr     = trapz(lc(:,1),lc(:,1).*lc(:,2));
087 scr     = trapz(lc(:,1),lc(:,1).^2.*lc(:,2));
088 scr     = sqrt(scr-mcr^2);
089 g       = lc2tr(lc,mcr,scr);
090 
091 f = [u u];
092 f(:,2) = tranproc(Z,fliplr(g));
093 
094 process = tranproc(L,f);
095 process = [(1:length(process))' process];
096 
097 %Check the result:
098 %save load.dat process -ascii
099 %[RFC TP]=getrfc;
100 %load rfc.dat
101 %param=[min(process(:,2)) max(process(:,2)) 100];
102 %fr=cc2fr(param,RFC);
103 %NT=fr2nt(fr);
104 %mu=pickdiag(NT);
105 
106 %Check the result without reference to getrfc:
107 LCe = dat2lc(process);
108 max(lc(:,2))
109 max(LCe(:,2))
110 
111 clf
112 plot(lc(:,1),lc(:,2)/max(lc(:,2)))
113 hold on
114 plot(LCe(:,1),LCe(:,2)/max(LCe(:,2)),'-.')
115 title('Relative crossing intensity')
116 
117 %% Plot made by the function funplot_4, JE 970707
118 %param = [min(process(:,2)) max(process(:,2)) 100];
119 %plot(lc(:,1),lc(:,2)/max(lc(:,2)))
120 %hold on
121 %plot(levels(param),mu/max(mu),'--')
122 %hold off
123 %title('Crossing intensity')  
124 %watstamp;
125 
126 % Temporarily
127 %funplot_4(lc,param,mu);
128

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