%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Ricardo Carretero 08/Dec/2003
% Integrates a SYSTEM of 4 PDEs eqs in 1D using finite differences
% modelling the behaviour of bridges:
%
%  theta_tt - e1*theta_xx = -(6*K/m)*cos(theta)*sin(theta) -d1*theta_t + h1(x,t)
%  y_tt - e2*y_xxxx = -(2*K/m)*y +g-d2*y_t + h2(x,t)
%
% split into: (u=y-gm/(2K))
%
%  theta_t = omega
%  u_t = v
%  omega_t = e1*theta_xx -(6*K/m)*cos(theta)*sin(theta) -d1*omega + h1(x,t)
%  v_t = e2*u_xxxx -(2*K/m)*u -d2*v + h2(x,t)
%
% Ex:  bridge(N,dt,tf,disp,savesnaps)
%
% npts=N, dt:time step, tf:t_final, domain:[0,L]
%
%Ex: 
%global allx allt allu allv alltheta allomega;
%global e1 e2 d1 d2 L l K m
%global dxx sixKm twoKm u0 theta0 sintwopiLx sinpiLx;
%N=50;dt=0.001;tf=50;snaps=100;disp=50;
%bridge(N,dt,tf,disp,snaps)
%save filename
%load filename
%bridgefilm
%

function [] = bridge(N,dt,tf,exposh,saveframes)

 global allx allt allu allv alltheta allomega;
 global e1 e2 d1 d2 L l K m
 global dxx sixKm twoKm u0 theta0 sintwopiLx sinpiLx;
 
 ran=1;
 noise=0;
 L=1000;
 d1=0.01;
 d2=0.01;
 e1=0.01;
 e2=0.0001;
 K=1000;
 m=2500;
 l=6;
 
 twoKm=2*K/m;
 sixKm=6*K/m;
 
 format compact
 maxstep=tf/dt;
 dispscreen= fix(maxstep/10); 

 if(exposh>0) stopdisp= fix(maxstep/exposh); 
 else stopdisp=9e+99;
 end
 if(saveframes>0)stopsave= fix(maxstep/saveframes);
 else stopsave=9e+99;
 end
 
 % Spatial computational domain defined

 dx=L/(N-1);
 x=linspace(0,L,N);
 dxx = 1/dx^2;
 sinpiLx=sin(pi*x/L);
 sintwopiLx=sin(2*pi*x/L);
 
 rand('state',0)       %this fixes the seed for the random generation
 
 strlength=0;
 fs=16;

 allt = 0:stopsave*dt:tf;
 allu = zeros(saveframes+1,N);
 allv = zeros(saveframes+1,N);
 alltheta = zeros(saveframes+1,N);
 allomega = zeros(saveframes+1,N);
 allx = x;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 fprintf('\nL \t\t= %f\n',L);
 fprintf('l \t\t= %f\n',l);
 fprintf('dt \t\t= %f\n',dt);
 fprintf('dx \t\t= %f\n',dx);
 fprintf('N \t\t= %f\n',N);
 fprintf('d1 \t\t= %f\n',d1);
 fprintf('d2 \t\t= %f\n',d2);
 fprintf('e1 \t\t= %f\n',e1);
 fprintf('e2 \t\t= %f\n',e2);
 fprintf('m \t\t= %f\n',m);
 fprintf('K \t\t= %f\n',K);
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% GENERATING ICs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 u0 = 0;
 theta0 = 0; 
 
 nn=[1 4 3];		% which modes I am going to intitiate in u
 mm=[2 4 3];		% which modes I am going to intitiate in theta
 An=[1   0*0.1 0];		% the initial amplitudes for each of the u modes
 Am=[0.1 0*0.02 0];		% the initial amplitudes for each of the theta modes
 
 vini = zeros(1,N);
 omegaini = zeros(1,N);
 
 uini=u0;
 thetaini=theta0;
 
 for ii=1:length(nn)
  uini = uini + An(ii)*sin(nn(ii)*pi*x/L); 
 end
 for ii=1:length(mm)
  thetaini = thetaini + Am(ii)*sin(mm(ii)*pi*x/L); 
 end
 
 % add noise:
 uini    =uini    +0*0.001*(rand(1,N)-0.5);
 thetaini=thetaini+0*0.001*(rand(1,N)-0.5);
 
 save=0;
 t=0;
 u=uini;
 v=vini;
 theta=thetaini;
 omega=omegaini;

 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot IC
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 bridgeplot(x,u,v,theta,omega,N,1); 
 axis tight;
 ax=axis;
 MM=max(abs(ax(5)),abs(ax(6)));
 ax(5)=-MM;
 ax(6)= MM;
 axis(ax);
 mytext=['T=(0)/(' num2str(tf) ')'];
 pos=[L*0.8,l*0.8,MM];
 text(pos(1),pos(2),pos(3),mytext);
 pause(0.0001)
 fprintf('Press any key\n')
 pause
 
 save=save+1;
 allu(save,:) = u;
 allv(save,:) = v;
 alltheta(save,:) = theta;
 allomega(save,:) = omega;
     
 fprintf('\nt=%5.2f ',t);
 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Main loop
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 for kk=1:maxstep
    
  t=t+dt;
  uN = u+dt*deriv_u(u,v,theta,omega,x,t,N);
  thetaN = theta+dt*deriv_theta(u,v,theta,omega,x,t,N);
  vN = v+dt*deriv_v(u,v,theta,omega,x,t,N);
  omegaN = omega+dt*deriv_omega(u,v,theta,omega,x,t,N);

  %Fixed BCs
  uN(1)=u0;
  uN(N)=u0;
  vN(1)=0;
  vN(N)=0;
  thetaN(1)=theta0;
  thetaN(N)=theta0;

  u = uN;
  v = vN;
  theta = thetaN;
  omega = omegaN;
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% include noise
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   if(noise>0)
    u  = u+noise*(rand(N)-0.5);
   end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display time on screen
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   if (fix(kk/dispscreen)==kk/dispscreen)
    fprintf('%5.2f ',t);
    strlength=strlength+1;
    if strlength>=10
     strlength=0;
     fprintf('\n');
    end
   end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% store solution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   if (fix(kk/stopsave)==kk/stopsave)
    save=save+1;
    allu(save,:) = u;
    allv(save,:) = v;
    alltheta(save,:) = theta;
    allomega(save,:) = omega;
   end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   if(fix(kk/stopdisp)==kk/stopdisp&kk>1)
     bridgeplot(x,u,v,theta,omega,N,1); 
     mytext=['T=(' num2str(kk*dt) ')/(' num2str(tf) ')'];
     axis(ax)
     text(pos(1),pos(2),pos(3),mytext)
     pause(0.0001)
 %     pause
   end

 end
 
 figure(3)
 subplot(2,1,1)
 set(gca,'FontSize',[fs]);
 surf(allx,allt,allu);
 shading interp
 xlabel('x'); ylabel('t'); zlabel('u(x,t)');
 subplot(2,1,2)
 set(gca,'FontSize',[fs]);
 surf(allx,allt,alltheta);
 shading interp
 xlabel('x'); ylabel('t'); zlabel('\theta(x,t)');
 pause(0.0001)

 fprintf('\n');
 return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function du = deriv_u(u,v,theta,omega,x,t,N);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 global e1 e2 d1 d2 L l K m
 global dxx sixKm twoKm u0 theta0 sintwopiLx sinpiLx;

% Fixed BCs
 v(1)=0;
 v(N)=0;
 du = v;

 return;
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dv = deriv_v(u,v,theta,omega,x,t,N);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 global e1 e2 d1 d2 L l K m
 global dxx sixKm twoKm u0 theta0 sintwopiLx sinpiLx;

% Fixed BCs
 up = [ u0 u(1:N-1) ];
 um = [ u(2:N) u0 ];
 uxx=dxx*(up-2*u+um);
 uxxp = [ uxx(1) uxx(1:N-1) ];
 uxxm = [ uxx(2:N) uxx(N) ];
 uxxxx=dxx*(uxxp-2*uxx+uxxm);
 
 dv = -e2*uxxxx -twoKm*u -d2*v + h2(x,t,N);

 return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dtheta = deriv_theta(u,v,theta,omega,x,t,N);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 global e1 e2 d1 d2 L l K m
 global dxx sixKm twoKm u0 theta0 sintwopiLx sinpiLx;

% Fixed BCs
 omega(1)=0;
 omega(N)=0;
 dtheta = omega;

 return;
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function domega = deriv_omega(u,v,theta,omega,x,t,N);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 global e1 e2 d1 d2 L l K m
 global dxx sixKm twoKm u0 theta0 sintwopiLx sinpiLx;

% Fixed BCs
 thetap = [ theta0 theta(1:N-1) ];
 thetam = [ theta(2:N) theta0 ];
 thetaxx=dxx*(thetap-2*theta+thetam);
 
 %thetaxx(1)=2*thetaxx(2)-thetaxx(3);
 %thetaxx(N)=2*thetaxx(N-1)-thetaxx(N-2);
 
 %domega = e1*thetaxx -sixKm*theta -d1*omega + h1(x,t);
 domega = e1*thetaxx -sixKm*cos(theta).*sin(theta) -d1*omega + h1(x,t,N);

 return;
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function hh1 = h1(x,t,N);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 global e1 e2 d1 d2 L l K m
 global dxx sixKm twoKm u0 theta0 sintwopiLx sinpiLx;

 %hh1 = 1*0.06*sin(1.2*t).*sintwopiLx;
 hh1 = 0*0.06*sin(1.2*t).*sinpiLx;

 return;
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function hh2 = h2(x,t,N);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 global e1 e2 d1 d2 L l K m
 global dxx sixKm twoKm u0 theta0 sintwopiLx sinpiLx;

 hh2 = 0;
 return;
