% The script file pwsemg04 computes examples of twice
% differentiable realizations with a dominant frequency.

Tp=1; omega=2*pi/Tp; dt=0.02; n=2400; eta=1/1; omegadt=omega*dt;
etadt=eta*dt; gamma=exp(-etadt); sigma2=sqrt(1); sigma0=sqrt(8/3);
alfa=sigma0*sqrt(2*eta);

% The fundamental solution

ME=eye(2); M0=zeros(2);
MR=[cos(omegadt),sin(omegadt);-sin(omegadt),cos(omegadt)];

FI1=gamma*[MR,M0,M0];
FI2=gamma*[etadt*MR,MR,M0];
FI3=gamma*[etadt*etadt/2*MR,etadt*MR,MR];
FI=[FI1;FI2;FI3]; 


% The asymptotic covariance matrix P

PP=[ME,ME/2,ME/4;ME/2,ME/4,ME*3/8;ME/4,ME*3/8,ME*3/8];
PP=sigma2*sigma2*sigma0*sigma0*PP;

% The triangular matrix in initial conditions

PA=[1,1/2,1/4;1/2,1/2,3/8;1/4,3/8,3/8];
p11=sqrt(PA(1,1));p12=0;p13=0;p21=PA(1,2)/p11;p31=PA(1,3)/p11; 
p22=sqrt(PA(2,2)-p21*p21);p23=0; p32=(PA(2,3)-p31*p21)/p22;
p33=sqrt(PA(3,3)-p31*p31-p32*p32);
p=[p11,p12,p13;p21,p22,p23;p31,p32,p33];
p=sigma2*sigma0*p;

p1=[p(1,1)*ME,p(1,2)*ME,p(1,3)*ME];
p2=[p(2,1)*ME,p(2,2)*ME,p(2,3)*ME];
p3=[p(3,1)*ME,p(3,2)*ME,p(3,3)*ME];
pp=[p1;p2;p3;];

% The triangular matrix in the simulation in one step

FIA=gamma*[1,0,0;etadt,1,0;etadt*etadt/2,etadt,1]; 
QQA=PA-FIA*PA*FIA';
qa11=sqrt(QQA(1,1)); qa21=QQA(1,2)/qa11; qa31=QQA(1,3)/qa11; 
qa12=0; qa13=0;qa22=sqrt(QQA(2,2)-qa21*qa21);qa23=0; 
qa32=(QQA(2,3)-qa31*qa21)/qa22;
qa33=sqrt(QQA(3,3)-qa31*qa31-qa32*qa32);
qa=[qa11,qa12,qa13;qa21,qa22,qa23;qa31,qa32,qa33]; 

qq=[qa11*ME,qa12*ME,qa13*ME;qa21*ME,qa22*ME,qa23*ME];
qq=[qq;qa31*ME,qa32*ME,qa33*ME];

U=[randn,randn,randn,randn,randn,randn];
ZD0=pp*U'; MZD=[ZD0'];

for r=1:n
   
   U=[randn,randn,randn,randn,randn,randn];
   ZD1=FI*ZD0+qq*U';
   MZD=[MZD;ZD1'];
   ZD0=ZD1;
end   

T=1:n+1; T=T*dt;
X=MZD(:,5); Y=MZD(:,6); AZ=abs(X+i*Y);

plot(T,X,'r',T,AZ,'k',T,-AZ,'k')
pause
close
   
U=zeros(1,6);
ZD0=pp*U'; MZD=[ZD0'];

for r=1:n
   
   U=[randn,randn,randn,randn,randn,randn];
   ZD1=FI*ZD0+qq*U';
   MZD=[MZD;ZD1'];
   ZD0=ZD1;
end   

T=1:n+1; T=T*dt;
X=MZD(:,5); Y=MZD(:,6); AZ=abs(X+i*Y);

plot(T,X,'r',T,AZ,'k',T,-AZ,'k')
pause
close
   
