% notopo.m
% Mousumi Roy - 10/15/02
% All quantities are in SI units

% fixed constants, including acceleration due to gravity
rhoc = 2600;
rhom = 3300;
E  = 70e9;
nu = 0.25;
Te = 20e3;
g  = 10;

% find delrho, the density contrast for crust and mantle
delrho = (rhom - rhoc);

% find flexural rigidity
D = (E*Te^3)/(12*(1-nu^2));

% find alpha, the flexural parameter
alpha = [4*D/(delrho*g)]^(0.25);

% User input here the appropriate loads
V0 = input('What is end load Vo ?  ');
M0 = input('What is bending moment Mo ?  ');
col = input('What line color (r,g,b,m,c,y)?   ', 's');

% with no topo loading, the constants in the solution to 
% the flexure equation, A and B, are:

A = (V0*alpha+M0)*(alpha^2)/(2*D);
B = (-M0)*(alpha^2)/(2*D);

% now find the deflection for position x from 0 to 800 km
x = [0:10000:800e3];

w = exp(-x./alpha).*(B*sin(x./alpha) + A*cos(x./alpha));

% now output graphics

figure(1)
     plot(x/1000, w/1000, 'color',col);
     xlabel('Position, km')
     ylabel('Deflection, km')

     disp('Note: to continue plotting in the same axes, type: hold on')
     disp('To save the figure as filename.ps, type: print filename.ps -dpsc')
     disp('Or, if you want to clear the figure and start again, type: clf')
     disp('To save the deflection curve into a new variable, w1, type: w1=w')
