(* ----------------------------------------------------------------------------------------- *)
(* This notebook solves the one dimensional wave propagation                                 *)
(* in a finite string. We used the seperation of variables technique.                        *)
(* The resulting sperate ODE in space and in time were solved seperately.                    *) 
(* Two Boundary conditions on fixed displacements and one on zero initial                    *) 
(* velocity resulted in a cos(λt) dependency on time and Sin(λx/c1) and sin(λ(x-α)/c2)       *)
(* in space. Finally the initial string displacements was used to determine the coefficients *)
(* ----------------------------------------------------------------------------------------- *)

(* --------------------------------- *)
(*                Data               *)
(* --------------------------------- *)
     (Subscript[c, 1])=20;     (Subscript[c, 2])=40;      L=1;      α=0.5*L;                  

(* ------------------------------------- *)
 (* Root finding Range and Step Size Data *)
(* ------------------------------------- *)
                min = 0.1;max =1000; inc = 1.0;

(* ------------------------------------------------------------------------------ *)
(* Nonlinear function to determine the eigenvalues λ based on Boundary Conditions *)
(* ------------------------------------------------------------------------------ *)
                   F[λ_]:=(1/(Subscript[c, 2]))*Sin[((λ *α)/(Subscript[c, 1]))]*Cos[((λ *(L-α))/(Subscript[c, 2]))]+(1/(Subscript[c, 1]))* Cos[((λ *α)/(Subscript[c, 1]))]*Sin[((λ *(L-α))/(Subscript[c, 2]))];

(* --------------------------------------------------------------------------------*) 
(* Displacement Continuity at x = α                                                *)
(* ------------------------------------------------------------------------------- *)
us[x_,λ_,α_,L_]:= Piecewise[{{Sin[((λ *(L-α))/(Subscript[c, 2]))]*Sin[(( λ *x)/(Subscript[c, 1]))],0≤x≤α},{Sin[((λ *α)/(Subscript[c, 1]))]*Sin[((λ *(L-x))/(Subscript[c, 2]))], α≤x≤L} } ];

(* ------------------------------------------------------------------------------- *)
(* Slope Continuity at x = α                                                       *)
(* ------------------------------------------------------------------------------- *)
uc[x_,λ_,α_,L_]:= Piecewise[{{-(1/(Subscript[c, 2]))*Cos[((λ *(L-α))/(Subscript[c, 2]))]*Sin[((λ *x)/(Subscript[c, 1]))],0≤x≤α},{(1/( (Subscript[c, 1])))*Cos[((λ *α)/(Subscript[c, 1]))]*Sin[((λ *(L-x))/(Subscript[c, 2]))],α≤x≤L}}];

(* ----------------------------------------------------------------------- *)
(*                Displacement Space-Time Solution                     *)
(* ----------------------------------------------------------------------- *)
u[x_,t_,λ_,α_,L_]:= Piecewise[{{us[x,λ,α,L]*Cos[λ*t],Abs[Sin[((λ *(L-α))/(Subscript[c, 2]))]]≥Sqrt[$MachineEpsilon] },
                                                                      {uc[x,λ,α,L]*Cos[λ*t], Abs[Sin[((λ *(L-α))/(Subscript[c, 2]))]]<Sqrt[$MachineEpsilon] }}];

(* --------------------------------------------------------------------------- *)
(*                           Initial Diaplacement                              *)   
(* --------------------------------------------------------------------------- *)                                          
g[x_,α_,L_]:=Sin[π*x];

(* ---------------------------------------------- *)
(*               Weight Function                  *)   
(* ---------------------------------------------- *)
w[c1_,c2_,α_]:=Piecewise[  {  {(1/c1)*(1/c1),0≤ x≤ α},{(1/c2)*(1/c2),α≤ x≤ L}  }  ];

(* ------------------------------------------------------------------ *)
(* Find all eigen values between min and max using inc as a step size *)
(* ------------------------------------------------------------------ *)
λn={};
For[ i=min,i < max, i+=inc, 
        If[F[i]*F[i+inc]≤0,roots=FindRoot[F[λ],{λ,i}]; λn=Append[λn,roots[[1,2]]],False]];

(* ----------------------------------------- *)
(* Print and plot out all the eigenvalues λn *)
(* ----------------------------------------- *)
λn
Plot[{ us[x,λn[[1]],α,L],us[x,λn[[2]],α,L],us[x,λn[[3]],α,L],
        us[x,λn[[4]],α,L],us[x,λn[[5]],α,L],us[x,λn[[6]],α,L]  },
           {x,0,L},Filling->Axis];

(* --------------------------------------- *)
(* Plot Initial Displacement of the string *)
(* --------------------------------------- *)
Plot[g[x,α,L],{x,0,L},Filling->Axis];

(* ======================================================== *)
(*               (∑_(p=1)^N)An u(x,λn,t=0)   =             g(x)    *)
(*   u(x,λm,t=0) (∑_(p=1)^N)An u(x,λn,t=0)   =  u(x,λm,t=0)g(x)    *)
(* (∫_0^L)u(x,λm,t=0) (∑_(p=1)^N)An u(x,λn,t=0)\[DifferentialD]x = (∫_0^L)u(x,λm,t=0)g(x)\[DifferentialD]x *)
(*        Amatrix . Unknown An       =        RHS b        *)
(* ======================================================== *)

(* ------------------------------------------------ *)
(* Compute, print, LogPlot Diagonal Amatrix         *)
(* ------------------------------------------------ *)
A=Range[Length[λn]];
For[i=1,i≤ Length[λn],i++,
        A[[i]]=(∫_0^L)u[x,0,λn[[i]],α,L]*u[x,0,λn[[i]],α,L]*w[(Subscript[c, 1]),(Subscript[c, 2]),α]\[DifferentialD]x ];
A
ListLogPlot[A,PlotRange-> {10^(-10),10^(-2)},Filling->Axis]

(* ----------------------------------------- *)
(* Compute, print, LogPlot RHS b             *)
(* ----------------------------------------- *)
b=Range[Length[λn]];
For[i=1,i≤ Length[λn],i++,
        b[[i]]=(∫_0^L)u[x,0,λn[[i]],α,L]*g[x,α,L]*w[(Subscript[c, 1]),(Subscript[c, 2]),α]\[DifferentialD]x ];
b
ListLogPlot[A,PlotRange-> {10^(-10),10^(-2)},Filling->Axis];

(* ----------------------------------------------------- *)
(* Compute Fourier coefficients An corresponding to g(x) *)
(* ----------------------------------------------------- *)
An=b/A;
ListLogPlot[Abs[An],PlotRange-> {10^(-32),10},Filling->Axis];

(* ============================================== *)
(* Construct the space-time displacement solution *)
(*   (∑_(p=1)^N)An u(x,λn,t)                              *)
(* ============================================== *)
U[x_,t_]:=(∑_(p=1)^(Length[λn]))An[[p]]*u[x,t,λn[[p]],α,L]
(* ----------------------------------------------------- *)
(* Plot Fourier approximation of Initial conditions g(x) *)
(* ----------------------------------------------------- *)
Plot[U[x,0],{x,0,L},Filling->Axis,PlotRange-> All ]

(* --------------------------------------------------------------------------------- *)
(*                     Visualize space time solution in 2D                           *)
(* --------------------------------------------------------------------------------- *)
Manipulate[Plot[{g[x,α,L],U[x,T]},{x,0,L},Filling->Axis,PlotRange-> {-1.1,1.1},Frame->True ],{T,0.000,1.0,0.001}]

(* ----------------------------------- *)
(* Visualize space time solution in 3D *)
(* ----------------------------------- *)
Plot3D[U[x,T],{x,0,L},{T,0,0.1},
              PlotRange -> {-1.5,1.5},Background->Black,Mesh->None,PlotPoints->40,AxesStyle->White,Boxed->False ]

Piecewise[  {        {0,0≤ x≤ (α/4)},{1,(α/4)≤ x≤ (α/2)},{0,(α/2)≤ x≤ L}       }  ]

Plot[Sin[π*x],{x,0,1}]

α

Sin[10(π/10+-α/3)]

Solve[10(x-α/3)==π]

Sin[x/π]