Home > wafo > wdemos > itmkurs > itmkurs_lab3.m

itmkurs_lab3

PURPOSE ^

Script to computer exercises 3

SYNOPSIS ^

This is a script file.

DESCRIPTION ^

 ITMKURS_LAB3 Script to computer exercises 3
  
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  Power Spectrum and Rainflow Analysis
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 %ITMKURS_LAB3 Script to computer exercises 3
002 % 
003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
004 % Power Spectrum and Rainflow Analysis
005 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
006 
007 %
008 %%Power Spectrum and Rainflow Analysis
009 %%Power Spectrum
010 %
011 
012 S=jonswap;
013 S=oscspec;
014 S=torsethaugen;
015 load deep.dat, S=dat2spec(deep); 
016 
017 %Solution Problem1:
018 
019 L=spec2mom(S,4);
020 f0=sqrt(L(2)/L(1))/2/pi
021 s=sqrt(L(1)) 
022 alfa=f0/(sqrt(L(3)/L(2))/2/pi)
023 wspecplot(S);
024 2*pi*f0
025 
026 beta=3:0.1:5;
027 gam=5E-9;
028 dpl=gam*spec2dplus(S,beta);
029 Tfpl=(1./dpl)/3600/24  %in days
030 clf, plot(beta,Tfpl)
031 
032 %
033 %%Simulation of $X$ and the damage intensity
034 %
035 
036 
037 T=1000/f0 % in seconds
038 gam=5E-9;
039 100*T*gam*spec2dplus(S,4.22) % in percent
040 
041 
042 max(S.w)/pi
043 dt=1./f0/60
044 S1=specinterp(S,dt);
045 [max(S1.w) pi/dt]
046 clf, wspecplot(S), hold on, wspecplot(S1,1,'r.')
047 clf, plot(beta,Tfpl)
048 hold on
049 for  i=1:5
050    X=spec2sdat(S1,60000);
051    tp=dat2tp(X);
052    rfc=tp2rfc(tp);
053    db=cc2dam(rfc,beta);
054    plot(beta,(1000/f0)*(1./db)*(1/gam)/3600/24,'r.') % Tf in days
055 end
056 
057 %
058 %%Theoretical computation of damage intensity
059 %
060 
061 
062 S=createspec;
063 w=levels([0 4 253]);
064 S.S=2*exp(-2*abs(w-2).^1.4);
065 S.S=S.S+9*exp(-8*(w-0.5).^2);
066 S.w=w; 
067 S.note=['S(w)=9*exp(-8(w-0.5)^2)+2*exp(-2|w-2|.^1.4)']
068 wspecplot(S)
069 
070 L=spec2mom(S,4);
071 f0=sqrt(L(2)/L(1))/2/pi
072 dpl=gam*spec2dplus(S,beta);
073 Tfpl=(1./dpl)/3600/24  %in days
074 dt=1./f0/60
075 S1=specinterp(S,dt);
076 figure(1)
077 clf, plot(beta,Tfpl)
078 hold on
079 for  i=1:10
080   X=spec2sdat(S1,60000);
081   tp=dat2tp(X);
082   rfc=tp2rfc(tp);
083   db=cc2dam(rfc,beta);
084   plot(beta,(1000/f0)*(1./db)*(1/gam)/3600/24,'r.') % Tf in days
085 end
086 
087 
088 
089 a=sqrt(L(3)/L(2))/2/pi
090 s=sqrt(L(1))
091 help spec2cmat
092 paramu=[-4.5*s 4.5*s 46];
093 nit=2;
094 figure(2), clf
095 [frfc fMm]=spec2cmat(S,[],'rfc',[],paramu,nit);
096 hold on, plot(rfc(:,2),rfc(:,1),'.')
097 clf, pdfplot(fMm)
098 dg=a*cmat2dam(paramu,frfc.f,beta);
099 figure(1), plot(beta,(1./dg)*(1/gam)/3600/24,'g') % Tf in days
100 
101 
102 
103 figure(2), clf
104 [frfc1 fMm1]=spec2cmat(S,[],'rfc',[],paramu,3);
105 dg1=a*cmat2dam(paramu,frfc1.f,beta);
106 figure(1)
107 plot(beta,(1./dg1)*(1/gam)/3600/24,'k') % Tf in days
108 
109 
110 help mctpsim
111 F{1,1}=fMm.f;
112 F{1,2}=fMm.f';
113 mctp=mctpsim(F,1000);
114 clf,plot(mctp(1:40)) 
115 
116 frfc1=zeros(size(frfc.f));
117 frfc1=triu(frfc.f,5);
118 fMm1=iter(frfc1,fMm.f,20,0.001);
119 clf,contour(fMm.f,20) 
120 hold on,  contour(fMm1,20)
121 F{1,1}=fMm1;
122 F{1,2}=fMm1';
123 mctp1=mctpsim(F,1000);
124 clf, subplot(2,1,1)
125 plot(mctp(1:40))
126 subplot(2,1,2)
127 plot(mctp1(1:40))
128 
129 
130 
131 
132 
133 
134 
135 
136

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