/*
	
	CantileverBeam Calculation Routine 
	
	* Methods  
		Concentrated(double,double)
		PartialUniform(double w,double a,double b)
		PartialInc(double w,double a,double b)
		PartialDec(double w,double a,double b)
		Couple(double mo, int ia)


*/

class CantileverBeam
{
	int noPoints=801;
	double E;			// modulus of elasiticiy
	double I;			// moment of intertia
	double L;			// span length
	int numPoints;		// number of points at which to compute values
	int i;
	// arrays to store distance, shear, moment and deflection
	double x[] = new double[noPoints];				// distance from the left end of the beam
	double v[] = new double[noPoints];				// shear force
	double m[] = new double[noPoints];				// moment
	double y[] = new double[noPoints];				// deflection
	double rr;			// reaction at right (fixed) support
	double mr;			// moment at right (fixed) supprot

	// CantileverBeam Constructor  
	// (E, I, Length of the beam, Number of points)
	CantileverBeam(double E,double I,double L, int n)
	{
		this.E = E;
		this.I = I;
		this.L = L;
		this.numPoints=n;

		for (i=0;i<numPoints;++i)
			x[i]=(double)i*L/((double)numPoints-1.);

	}




	// Compute concentrated load at distance a from the right end
	// Formulas 2
	void Concentrated(double w, double a)
	{

		// compute reaction
		double b=L-a;

		rr = w;

		mr = w*a;
		

		// compute shear, moemnt and deflection till the point a

		for (i=0;i<numPoints;++i)
		{
			if (x[i]<b)
			{
				v[i] = 0.;
			
				m[i] = 0.;
			
				y[i] = -rr/(6.*E*I)*(-a*a*a+3.*a*a*L-3.*a*a*x[i]);

			}

		    else
			{
				v[i] = -rr;
			
				m[i] = -rr*(x[i]-b);
			
				y[i] = -rr/(6*E*I)*(Math.pow(x[i]-b,3)-3*a*a*(x[i]-b)+2*a*a*a);
			}
		}

	}


	// Partial Uniform Load from a to b,
	void ComputeLinear(double w1,double w2,double a,double b)
	{

		// compute reaction
    	rr = w1*(b-a);

	    mr = 0.5*rr*(a+b);
		
		   // compute shear, moemnt and deflection 

    	for (i=0;i<numPoints;++i)
	    {
		    
		    if (x[i]<(L-b)) 
    		{
	    		v[i] = 0;
    
    	    	m[i] = 0;			

		    	y[i] = -rr/(24*E*I)*(4*(a*a+a*b+b*b)*(L-x[i])-a*a*a-a*b*b-a*a*b-b*b*b);
			}
    		else 
	    	{
		    	if (x[i]<(L-a))
			   	{
			    	v[i] = -rr/(b-a)*(x[i]-L+b);
    
	    			m[i] = -0.5*rr/(b-a)*Math.pow((x[i]-L+b),2);
    
	    			y[i] = -rr/(24*E*I)* ( 6*(a+b)*Math.pow(L-x[i],2) -4*Math.pow(L-x[i],3)
		    				+Math.pow(L-x[i]-a,4)/(b-a));
				}
				else
				{
					v[i] = -rr;

					m[i] = -0.5*rr*(2*x[i]-2*L+a+b);

					y[i] = -rr/(12*E*I)*(3*(a+b)*Math.pow(L-x[i],2)-2*Math.pow(L-x[i],3));
				}
			}
		}
		
        if (w2-w1!=0.)
        {
            double W = 0.5*(w2-w1)*(b-a);
            
    		rr += W;

	    	mr += 0.5*W*(2*a+b);
		
		    // compute shear, moemnt and deflection 

    		for (i=0;i<numPoints;++i)
	    	{
		
		    	if (x[i]<(L-b)) 
			    {
				    v[i] += 0;

    				m[i] += 0;			

	    			y[i] += -W/(60*E*I)*((5*b*b+10*a*b+15*a*a)*(L-x[i])
		    			-b*b*b-2*b*b*a-3*b*a*a-4*a*a*a);
			    }
			    else 
			    {
				    if (x[i]<(L-a))
				    {
					    v[i] += -W*Math.pow(x[i]-L+b,2)/Math.pow(b-a,2);

					    m[i] += -W/3*(Math.pow(x[i]-L+b,3)/Math.pow(b-a,2));

    					y[i] += -W/(60*E*I)*((20*a+10*b)*Math.pow(L-x[i],2)-10*Math.pow(L-x[i],3)
	    					+5*Math.pow(L-x[i]-a,4)/(b-a)-Math.pow(L-x[i]-a,5)/Math.pow(b-a,2));
				    }
    				else
	    			{
		    			v[i] += -W;

			    		m[i] += -W/3*(-3*L+3*x[i]+b+2*a);

				    	y[i] += -W/(6*E*I)*((2*a+b)*Math.pow(L-x[i],2)-Math.pow(L-x[i],3));
				    }
			    }
			}
		}		
	}
    



	// Moment applied at distance a from the right end
	// Formulas 10
	void Couple(double mo, double a)
	{

		// compute reaction
		
		double b=L-a;

		rr = 0;
		mr = -mo;
		
		// compute shear, moemnt and deflection till the point a

		for (i=0;i<numPoints;++i)
		{
			if (x[i]<b)
			{
				v[i] = 0.;
			
				m[i] = 0.;
			
				y[i] = mo*a/(E*I)*(L-0.5*a-x[i]);

			}

			else
			{
				v[i] = 0;
			
				m[i] = mo;
			
				y[i] = 0.5*mo/(E*I)*(Math.pow(x[i]-L+a,2)-2*a*(x[i]-L+a)+a*a);
			}
		}
		

	}
			
}

/*
	
	EndSupports.java

  * Methods
		Concentrated(double W, double a)
		PartialUniform(double w, double b, double a)
		PartailInc(double w,double a,double b)
		Couple(double mo, double a)


*/



// Main Class for class Both ends supported beam calculation

class EndSupports
{
	int noPoints=801;
	double E;			// modulus of elasiticiy
	double I;			// moment of intertia
	double L;			// span length
	int numPoints;		// number of points at which to compute values
	int i;

	// arrays to store distance, shear, moment and deflection
	double x[]=new double[noPoints];		// distance from the left end of the beam
	double v[]=new double[noPoints];		// shear force
	double m[]=new double[noPoints];		// moment
	double y[]=new double[noPoints];		// deflection
	double rr1;			// reaction at left support
	double rr2;			// reaction at right support

	// EndSupports Class constructor

	EndSupports(double E,double I,double L,int n)
	{
		this.L = L;
		this.E = E;
		this.I = I;
		this.numPoints=n;

		for (i=0;i<numPoints;++i)
			x[i]=(double)i*L/((double)numPoints-1.);

	}


	// compute concentrated load distance a from the left end (a starts from 0)
	// formulas 12
	void Concentrated(double W, double a)
	{

		// compute reaction
		double b=L-a;
		
		rr1 = W*b/L;
		rr2 = W*a/L;
		
		// compute shear, moemnt and deflection

		for (i=0;i<numPoints;++i)
		{
		
			if (x[i]<a) 
			{
				v[i] = W*b/L;

				m[i] = W*b/L*x[i];			

				y[i] = -W*b*x[i]/(6*E*I*L)*( 2*L*(L-x[i])-b*b-(L-x[i])*(L-x[i]) );
			}
			else 
			{
				v[i] = -W*a/L;

				m[i] = W*a/L*(L-x[i]);

				y[i] = -W*a*(L-x[i])/(6*E*I*L)*( 2*L*b-b*b-Math.pow(L-x[i],2) );
			}
		}
	}



// compute partial uniform load situation 
	// Fomulas 14 and 16
	void ComputeLinear(double w1,double w2, double a, double b)
	{
		double d;
		double c;
		double W;

		d = L-0.5*b-0.5*a; // distance from the right end to the center of the load
		c = b-a;		// length of load
        W = w1*c;
		rr1 = W * d / L;        // reaction at the left end
		rr2 = W/L*(a+0.5*c);	    // reaction at the right end
		for (i=0;i<numPoints;++i)
		{
		
			if (x[i]<a) 
			{
				v[i] = rr1;

				m[i] = rr1*x[i];			

				y[i] = 1/(48*E*I)*(8*rr1*(Math.pow(x[i],3)-L*L*x[i])
					+W*x[i]*(8*d*d*d/L-2*b*c*c/L+c*c*c/L+2*c*c));
			}
			else 
			{
				if (x[i]<b)
				{
					v[i] = rr1-W*(x[i]-a)/c;

					m[i] = rr1*x[i]-W*Math.pow(x[i]-a,2)/(2*c);

					y[i] = 1/(48*E*I)*(8*rr1*(Math.pow(x[i],3)-L*L*x[i])
						+W*x[i]*(8*d*d*d/L-2*b*c*c/L+c*c*c/L+2*c*c)
						-2*W*Math.pow(x[i]-a,4)/c);
				}
				else
				{
					v[i] = rr1-W;

					m[i] = rr1*x[i]-W*(x[i]-0.5*a-0.5*b);

					y[i] = 1/(48*E*I)*(8*rr1*(Math.pow(x[i],3)-L*L*x[i])
						+W*x[i]*(8*d*d*d/L-2*b*c*c/L+c*c*c/L)
						-8*W*Math.pow(x[i]-0.5*a-0.5*b,3)+W*(2*b*c*c-c*c*c));
				}
			}
		}
	

		


	// Compute partial triangular load from a to b, load varying from zero to w
	// a and b measured from the left and. 
	// formulas 16
	

		// compute reaction
		
		if(w2-w1!=0)
		{

    		d = L-a/3.-2./3.*b; // distance from the right end to the center of the load
    		W = 0.5*(w2-w1)*c;		// total load

    		double rr1a = W*d/L;
    		rr1+= rr1a;
	    	rr2 += W*(L-d)/L;
		
		    // compute shear, moment and deflection 

    		for (i=0;i<numPoints;++i)
	    	{
		
		    	if (x[i]<a) 
			    {
				    v[i] += rr1a;

    				m[i] += rr1a*x[i];			

	    			y[i] += 1./(6.*E*I)*( rr1a*(Math.pow(x[i],3.)-L*L*x[i])
		    			    +W*x[i]*(d*d*d/L+1./6.*c*c*(1.-b/L)+17./270.*c*c*c/L));
			    }
    			else 
	    		{
		    		if (x[i]<b)
			    	{
				    	v[i] += rr1a-Math.pow((x[i]-a)/c,2)*W;

					    m[i] += rr1a*x[i]-W*Math.pow(x[i]-a,3)/(3*c*c);						

    					y[i] += 1./(6.*E*I)*(rr1a*(Math.pow(x[i],3.)-L*L*x[i])
	    					    -W/10.*Math.pow(x[i]-a,5.)/(c*c)
		    				    +W*x[i]*(d*d*d/L+c*c/6.-c*c*b/6./L+17./270.*c*c*c/L));
			    	}
				    else
				    {
					    v[i] += rr1a-W;

    					m[i] += rr1a*x[i]-W/3*(3*x[i]-a-2*b);
    
	    				y[i] += 1./(6.*E*I)*( rr1a*(Math.pow(x[i],3.)-L*L*x[i])
		    				    -W* ( Math.pow(x[i]-a/3.-2*b/3.,3)-d*d*d*x[i]/L
			    			    -b/6.*c*c*(1.-x[i]/L)+17./270.*c*c*c*(1-x[i]/L) ) );
				    }
			    }
		    }

		}

	}
	// Moment applied mo at distance a from the left end
	// formulas 20
	void Couple(double mo, double a)
	{

		// compute reaction
		rr1 = -mo/L;
		rr2 = -rr1;
		
		// compute shear, moemnt and deflection till the point a

		for (i=0;i<numPoints;++i)
		{
		
			if (x[i]<a) 
			{
				v[i] = rr1;

				m[i] = rr1*x[i];			

				y[i] = mo/(6*E*I)*( (6*a-3*a*a/L-2*L)*x[i]-Math.pow(x[i],3)/L );
			}
			else 
			{
				v[i] = rr1;

				m[i] = rr1*x[i]+mo;

				y[i] = mo/(6*E*I)*( 3*a*a+3*x[i]*x[i]-Math.pow(x[i],3)/L
					-(2*L+3*a*a/L)*x[i] );
			}
		}
		

	}

}


/*
	
	FixedEnds.java: Applet 

	* Methods
		Concentrated(double W, double a)
		PartialUniform(double w, double b, double a)
		PartialInc(double w,double a,double b)
		Couple(double mo, double a)


*/


class FixedEnds
{
	int noPoints=801;
	double E;			// modulus of elasiticiy
	double I;			// moment of intertia
	double L;			// span length
	int numPoints;		// number of points at which to compute values
	int i;

	// arrays to store distance, shear, moment and deflection
	double x[]=new double[noPoints];		// distance from the left end of the beam
	double v[]=new double[noPoints];		// shear force
	double m[]=new double[noPoints];		// moment
	double y[]=new double[noPoints];		// deflection
	double rr1;			// reaction at left fixed end
	double rr2;			// reaction at right fixed end
	double mm1;			// moment at left end
	double mm2;			// moment at right end
	
	// FixedEnds Class constructor

	FixedEnds(double E,double I,double L,int n)
	{
		this.L = L;
		this.E = E;
		this.I = I;
		this.numPoints=n;

		for (i=0;i<numPoints;++i)
			x[i]=i*L/(numPoints-1.);

	}


	// compute concentrated load wihch is at distance a from the left end
	// formulas 32
	void Concentrated(double W, double a)
	{
		// compute reaction
		double b=L-a;
		
		rr1 = W*b*b/(L*L*L)*(3*a+b);
		rr2 = W*a*a/(L*L*L)*(3*b+a);
		mm1 = W*a*b*b/(L*L);
		mm2 = W*a*a*b/(L*L);
		
		// compute shear, moemnt and deflection

		for (i=0;i<numPoints;++i)
		{
		
			if (x[i]<a) 
			{
				v[i] = rr1;

				m[i] = -W*a*b*b/L/L+rr1*x[i];			

				y[i] = W*b*b*x[i]*x[i]/(6*E*I*L*L*L)*(3*a*x[i]+b*x[i]-3*a*L);
			}
			else 
			{
				v[i] = rr1-W;

				m[i] = -W*a*b*b/L/L+rr1*x[i]-W*(x[i]-a);			

				y[i] = W*a*a*Math.pow(L-x[i],2)/(6*E*I*L*L*L)
					*( (3*b+a)*(L-x[i])-3*b*L );
			}
		}
	}
	// compute partial uniform load situation Fomulas 34
	void ComputeLinear(double w1,double w2, double a, double b)
	{
		double d;
		double c;
		double W;

		d = L-a/2.-b/2.; // distance from the right end to the center of the load
		c = b-a;		// length of load
		W = w1*c;		// total load

		rr1 = W/(4.*L*L)*(12.*d*d-8.*d*d*d/L+2.*b*c*c/L-c*c*c/L-c*c);        // reaction at the left end
		rr2 = W-rr1;
		mm1 = -W/24./L*(24.*d*d*d/L-6.*b*c*c/L+3.*c*c*c/L+4.*c*c-24.*d*d);
		mm2 = W/24./L*(24.*d*d*d/L-6.*b*c*c/L+3.*c*c*c/L+2.*c*c-48.*d*d+24.*d*L);	    // reaction at the right end
		for (i=0;i<numPoints;++i)
		{
		
			if (x[i]<a) 
			{
				v[i] = rr1;

				m[i] = -mm1+rr1*x[i];			

				y[i] = 1./(6.*E*I)*( rr1*Math.pow(x[i],3.)-3.*mm1*x[i]*x[i]);
			}
			else 
			{
				if (x[i]<b)
				{
					v[i] = rr1-W*(x[i]-a)/c;

					m[i] = -mm1+rr1*x[i]-0.5*W*Math.pow(x[i]-a,2.)/c;

					y[i] = 1./(6.*E*I)*( rr1*Math.pow(x[i],3.)-3.*mm1*x[i]*x[i]
						-W/4.*Math.pow(x[i]-a,4.)/c);
				}
				else
				{
					v[i] = rr1-W;

					m[i] = -mm1+rr1*x[i]-W*(x[i]-L+d);

					y[i] = 1./(6.*E*I)*(rr2*Math.pow(L-x[i],3.)-3.*mm2*Math.pow(L-x[i],2.));
				}
			}
		}
	
	
	// Compute partial triangular load from b to a, load varying from zero to w
	// a and b measured from the right and. 
	// formulas 36
	

        if (w2-w1!=0.)
        {

    		d = L-a/3.-2./3.*b; // distance from the right end to the center of the load
    		W = 0.5*(w2-w1)*c;		// total load

	    	double rr1a = W/L/L*(3.*d*d-c*c/6.+b*c*c/3./L-17.*c*c*c/135./L-2.*d*d*d/L);        // reaction at the left end
		    rr1 += rr1a;
		    rr2 += W-rr1a;

            
    		double mm1a = -W/L*(d*d*d/L+c*c/9.+51./810.*c*c*c/L-c*c*b/6./L-d*d);
    		mm1 += mm1a;
    		mm2 += W/L*(d*d*d/L+c*c/18.+51./810.*c*c*c/L-c*c*b/6./L-2.*d*d+d*L);	    // reaction at the right end

	    	for (i=0;i<numPoints;++i)
		    {
		
    			if (x[i]<a) 
	    		{
		    		v[i] += rr1a;

			    	m[i] += -mm1a+rr1a*x[i];

				    y[i] += 1./(6.*E*I)*(rr1a*Math.pow(x[i],3.)-3.*mm1a*x[i]*x[i]);
    			}
	    		else 
		    	{
			    	if (x[i]<b)
				    {
					    v[i] += rr1a-W*Math.pow((x[i]-a)/c,2);

    					m[i] += -mm1a+rr1a*x[i]-W/3.*Math.pow(x[i]-a,3.)/(c*c);

	    				y[i] += 1./(E*I)*(1./6.*rr1a*Math.pow(x[i],3.)-.5*mm1a*x[i]*x[i]
		    				    -1./60.*W*Math.pow(x[i]-a,5.)/c/c );
			    	}
				    else
    				{
	    				v[i] += rr1a-W;

		    			m[i] += -mm1a+rr1a*x[i]-W*(x[i]-L+d);

			    		y[i] += 1./(6.*E*I)*(rr1a*(Math.pow(x[i],3.)-3.*L*L*x[i]+2*L*L*L)
				    		    -3.*mm1a*Math.pow(L-x[i],2.)+W*(3.*d*d*x[i]+d*d*d
					    	    -3.*d*d*L-Math.pow(x[i]-L+d,3.)));
				    }
    			}
	    	}
	    }


	}



	// Moment applied mo at distance a from the left end
	// formulas 37
	void Couple(double mo, double a)
	{

		// compute reaction
		rr1 = -6.*mo/L/L/L*(a*L-a*a);
		rr2 = -rr1;

		mm1 = -mo/L/L*(4.*a*L-3.*a*a-L*L);
		mm2 = mo/L/L*(2.*L*a-.3*a*a);
		
		// compute shear, moemnt and deflection till the point a

		for (i=0;i<numPoints;++i)
		{
		
			if (x[i]<a) 
			{
				v[i] = rr1;

				m[i] = -mm1+rr1*x[i];			

				y[i] = -1./6./E/I*( 3.*mm1*x[i]*x[i]-rr1*Math.pow(x[i],3.) );
			}
			else 
			{
				v[i] = rr1;

				m[i] = -mm1+rr1*x[i]+mo;			

				y[i] = 1./6./E/I*( (mo-mm1)*(3.*x[i]*x[i]-6.*L*x[i]+3.*L*L)
					-rr1*(3.*L*L*x[i]-Math.pow(x[i],3.)-2.*L*L*L));
			}
		}

	}

}

/*
	
	OneFixOneSupport.java

*/


class OneFixOneSupport
{
	int noPoints=801;
	double E;			// modulus of elasiticiy
	double I;			// moment of intertia
	double L;			// span length
	int numPoints;		// number of points at which to compute values
	int i;
	String buffer;

	// arrays to store distance, shear, moment and deflection
	double x[]=new double[noPoints];		// distance from the left end of the beam
	double v[]=new double[noPoints];		// shear force
	double m[]=new double[noPoints];		// moment
	double y[]=new double[noPoints];		// deflection
	double rr1;			// reaction at left support
	double rr2;			// reaction at right support
	double mm2;			// reaction moment ar right fixed end

	//  Class constructor

	OneFixOneSupport(double E,double I,double L,int n)
	{
		this.L = L;
		this.E = E;
		this.I = I;
		this.numPoints=n;

		for (i=0;i<numPoints;++i)
			x[i]=i*L/(numPoints-1.);

	}

      // compute concentrated load distance b from the left end
        // formulas 22
        void Concentrated(double W, double a)
        {

                // compute reaction
                double b=L-a;
                
                rr1 = 0.5*W*(3.*a*a*L-a*a*a)/(L*L*L);
                rr2 = W-rr1;
                mm2 = 0.5*W*(a*a*a+2.*a*L*L-3.*a*a*L)/(L*L);
                
                // compute shear, moemnt and deflection

                for (i=0;i<numPoints;++i)
                {
                
                        if (x[i]<b) 
                        {
                                v[i] = rr1;

                                m[i] = rr1*x[i];                        

                                y[i] = 1/(6*E*I)*(rr1*(Math.pow(x[i],3)-3*L*L*x[i])+3*W*a*a*x[i]);
                        }
                        else 
                        {
                                v[i] = rr1-W;

                                m[i] = rr1*x[i]-W*(x[i]-L+a);

                                y[i] = 1/(6*E*I)*(rr1*(Math.pow(x[i],3)-3*L*L*x[i])
                                        +W*(3*a*a*x[i]-Math.pow(x[i]-b,3)));
                        }
                }

        }

	// compute partial uniform load situation Fomulas 24
	void ComputeLinear(double w1, double w2,double a, double b)
	{
		double d;
		double c;
		double W;

		d = L-b; // distance from the right end to the center of the load
		c = b-a;		// length of load
		W = w1*c;		// total load

		rr1 = W/(8.*L*L*L)*(4.*L*(a*a+a*b+b*b)-a*a*a-a*b*b-a*a*b-b*b*b);        // reaction at the left end
		rr2 = W-rr1;
		mm2 = -rr1*L+0.5*W*(a+b);	    // reaction at the right end

		for (i=0;i<numPoints;++i)
		{
		
			if (x[i]<L-b) 
			{
				v[i] = rr1;

				m[i] = rr1*x[i];			

				y[i] = 1./(E*I)*( rr1*(Math.pow(x[i],3.)/6.-.5*L*L*x[i])
					+W*x[i]*(0.5*a*a+.5*a*c+c*c/6.));
			}
			else 
			{
				if (x[i]<L-a)
				{
					v[i] = rr1-W*(x[i]-d)/c;

					m[i] = rr1*x[i]-W*Math.pow(x[i]-d,2.)/(2.*c);

					y[i] = 1./(E*I)*( rr1*(Math.pow(x[i],3.)/6.-.5*L*L*x[i])
					    +W*x[i]*(.5*a*a+.5*a*c+c*c/6.)
						-W*Math.pow(x[i]-d,4.)/(24.*c) );
				}
				else
				{
					v[i] = rr1-W;

					m[i] = rr1*x[i]-W*(x[i]-d-c/2.);

					y[i] = 1./(E*I)*( rr1*(Math.pow(x[i],3.)/6.-.5*L*L*x[i]+L*L*L/3.)
						+W*( Math.pow(a+0.5*c,3.)/6.-0.5*Math.pow(a+0.5*c,2.)*L
					    -1./6.*Math.pow(x[i]-d-0.5*c,3)+0.5*Math.pow(a+0.5*c,2.)*x[i]) );
				}
			}
		}

	// Compute partial triangular load from b to a, load varying from zero to w
	// a and b measured from the right and. 
	// formulas 26


        if (w2-w1!=0.)
        {
    		d = L-b; // distance from the right end to the center of the load
    		W = 0.5*(w2-w1)*c;		// total load

    		double rr1a = W/(20.*L*L*L)*((10.*a*b+15.*a*a+5.*b*b)*L-4.*a*a*a-2.*a*b*b-3.*a*a*b-b*b*b);        // reaction at the left end
            rr1 += rr1a;
	    	rr2 += W-rr1a;
    		mm2 += -rr1a*L+W/3.*(2.*a+b);	    // reaction at the right end
	    	for (i=0;i<numPoints;++i)
    		{
			    if (x[i]<L-b) 
			    {
				    v[i] += rr1a;

    				m[i] += rr1a*x[i];

	    			y[i] += 1./(E*I)*( rr1a*(Math.pow(x[i],3.)/6.-.5*L*L*x[i])
		        		   +W*x[i]*(0.5*a*a+a*c/3.+c*c/12.));
			    }
    			else 
	    		{
		    		if (x[i]<(L-a))
			    	{
				    	v[i] += rr1a-W*Math.pow((x[i]-d)/c,2.);

    					m[i] += rr1a*x[i]-W/3.*Math.pow(x[i]-d,3.)/(c*c);

    					y[i] += 1./(E*I)*( rr1a*(Math.pow(x[i],3.)/6.-.5*L*L*x[i])
	    	    			    +W*x[i]*(0.5*a*a+a*c/3.+c*c/12.)
		    	    			-W*Math.pow(x[i]-d,5.)/(60.*c*c) );
    				}
				    else
    				{
	    				v[i] += rr1a-W;

    					m[i] += rr1a*x[i]-W*(x[i]-d-2.*c/3.);

    					y[i] += 1./(E*I)*( rr1a*(Math.pow(x[i],3.)/6.-.5*L*L*x[i]+L*L*L/3.)
	    				   	    +W*( 0.5*Math.pow(a+c/3.,2.)*x[i]-1./6.*Math.pow(x[i]-d-2./3.*c,3.)
		    				    +1./6.*Math.pow(a+1./3.*c,3.)-0.5*Math.pow(a+c/3.,2.)*L));
			    	}
			    }
			}
		}

	}


	// Moment applied mo at distance a from the left end
	// formulas 30
	void Couple(double mo, double a)
	{

		// compute reaction
		rr1 = -3./2.*mo/L*(L*L-a*a)/(L*L);
		rr2 = 3./2.*mo/L*((L*L)-a*a)/(L*L);
		mm2 = 0.5*mo*(1.-3.*a*a/(L*L));
		
		// compute shear, moemnt and deflection till the point a

		for (i=0;i<numPoints;++i)
		{
		
			if (x[i]<a) 
			{
				v[i] = rr1;

				m[i] = rr1*x[i];			

				y[i] = mo/(E*I)*( (L*L-a*a)/(4.*L*L*L)*(3.*L*L*x[i]-Math.pow(x[i],3.))-(L-a)*x[i]);
			}
			else 
			{
				v[i] = rr1;

				m[i] = rr1*x[i]+mo;

				y[i] = mo/(E*I)*( (L*L-a*a)/(4.*L*L*L)*(3.*L*L*x[i]-Math.pow(x[i],3.))-L*x[i]
					+0.5*(Math.pow(x[i],2.)+a*a));
			}
		}

	}
		

}