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

% fixed constants
rhoc = 2600; % crustal density
rhom = 3300; %mantle density
E  = 70e9;   % Youngs modulus
nu = 0.25;   % Poissons ratio
Te = 20e3;   % effective elastic plate thickness
g  = 10;     % accln due to gravity

% 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 topo
T = input('What is plateau height in meters?  ');
a = input('What is plateau width in meters?  ');
col = input('What line color (r,g,b,m,c,y)?   ', 's');

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

A = ;
B = ; 

% now find the deflection for position x from 0 to 800 km
x = [0:10000:800e3];
dummy = ones(size(x));
xin  = (x<=a).*x;
xin  = [ 0 xin(find(xin))];
xout = (x>a).*x;
xout = [ xout(find(xout))];

w0    = exp(-x./alpha).*(B*sin(x./alpha) + A*cos(x./alpha));
win   = (rhoc*T*0.5/delrho)*(2 - cos((xin+a)/alpha).*exp(-(a+xin)/alpha)...
       - cos((xin-a)/alpha).*exp((xin-a)/alpha));
wout  = (rhoc*T*0.5/delrho)*(cos((a-xout)/alpha).*exp((a-xout)/alpha) ...
       - cos((xout+a)/alpha).*exp(-(xout+a)/alpha));

wp    = [win wout];
wtot  = w0 + wp;
% choose deflection to be positive downwards, topo to be positive upward
wtot  = -wtot;
w     = wtot;
box = (x<=a).*dummy*T + (x>a).*dummy*0;

% now output graphics

     plot(x/1000, box/1000); hold on;
     plot(x/1000, w/1000, 'color',col); hold on
     xlabel('Position, km')
     ylabel('Deflection, km')

     disp('BEFORE RE-RUNNING THE PROGRAM, NOTE:')
     disp('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')



     figure(1)
