// Submitted by Dr. Triranjita Srivastava. Assistant Professor, Physics Dept., Kalindi College, University of Delhi

//Operating system: Windows 8
//SCILAB Ver: 5.5.2
//Objectiv: Lagrangian formulation in classical mechanics with constraints
//Example: Simple Pendulum of length L (m)operating in gravitational field . After applying Lagrangian formulation this problem reduces to a simple second order differential equation [(d^2 theta/dt^2)+(g/L)sin(theta)]=0. Here theta is angular displacement.
// We implemented ordinary differential equation (ODE) Solver to solve the second order differential equation
//We present plot of solution of angular displacement for t=0 to t=10 seconds 

clear 
clc
L=input ('Enter the length of pendulum (m) L = ')
g = 9.8                    //acceleration due to gravity (m/s^2)
k=g/L
theta=input('Enter the initial angular displacement (radian) at (t = 0) = ') ;        // Initial angular displacement at t = 0
dt=input('Enter initial d_theta/dt (radian) at (t = 0) = ') ;        // Initial boundary condition d_theta/dt at t = 0
////************************************************************************
//// Function declaration for ODE
////************************************************************************
t=linspace(0,10,200)
function dx=f(t,x,k)
    dx(1)=x(2)
    dx(2)=-k*sin(x(1))
endfunction
////************************************************************************
//// Solving second order differential equation by ODE  solver
////************************************************************************
    y=ode([theta;dt],0,t,f)
    ysol=y(1,:)
    ydotsol = y(2,:)

////************************************************************************
//// Plotting the solution (angular displacement (theta)and d_theta/dt)
////************************************************************************
scf()
title('Solution of Simple Pendulum', 'fontsize',5)
ylabel('Solution --->','fontsize',5)
xlabel('t (sec) ---> ','fontsize',5)
a=get("current_axes")       //get the handle of the newly created axes
a.font_size=4
t=get("hdl")                //get the handle of the newly created object
t.font_size=5
plot(t,ysol,'r','linewidth',3)
plot(t,ydotsol,'k','linewidth',3)
h1 = legend(['$\theta$'; '$d\theta/dt$'])
h1.font_size=3


