% Flex2D script for variable flexural rigitidy in x direction
% Centered Finite Difference Technique
% uses flexural equation from Turcotte and Schubert, 1982
% d2/dx2(D(x)(d2w/dx2))+ (pm-pc)gw = q(x)
% written by Jay Chapman, Nov. 2015

clear all
close all

%Grid size, number of nodes in x direction
sizex = 201;

%Node spacing (meters)
delta = 5000.0;

% profile distance in km
x = [0:delta:delta*(sizex-1)];

% Young Modulus in Pascals
E = 70.0e9;
% Poisson's ratio
v = 0.25;
% Gravity in m/s^2
g = 9.8;
% Density of the mantle kg/m^3
pm = 3300.0;
% Density of the material filling resultant depression kg/m^3
pfill = 2700;

% Loads heights (meters) and densities (kg/m^3)
hload = zeros(1,sizex); %initialize vector for load height
pload = zeros(1,sizex); %initialize vector for density of load
te = ones(1,sizex);  %initialize vector for elastic thickness

%Comment out loop below for constant Elastic thickness
te = 30000*te;


%********** This section is to assign elastic thickness **********
for i=1:sizex
    if i<70
        te(i)=15000;
    elseif i>=100 && i<=102
        te(i)=20000;
    elseif i>102
        te(i)=30000;
    end
end

%******** This section is to assign the load ****************
for i=1:sizex
    if i>40 && i<=50
        hload(i)=1000;
        pload(i)=2700;
    elseif i>50 && i<=75
        hload(i)=2000;
        pload(i)=2700;
    elseif i>75 && i<=85
        hload(i)=1000;
        pload(i)=2700;
    end
end

% Compute flexural rigidity
D = te.^3.0*E/(12.0*(1.0-v^2.0));

% Compute load
tload = (hload.*pload*g)';

% --------------------------------------------------------------------
% These are variables that go into making a sparse matrix for the inversion
% "sizeind" is the number of non-zero values in the sparse matrix
% "rowi" and "coli" point the position of the non-zero element "values"
sizeind = ((sizex-4)*5)+4; 
rowi = zeros(sizeind,1);
coli = zeros(sizeind,1);
values = zeros(sizeind,1);



% These are the outer non-zero terms
count = 1;
for i = [1 2 sizex-1 sizex]
       %sets the diagonal on first two and last two rows in DM =1
       %corresponds to the first and last two entries in the sizex vector
        rowi(count)=i;
        coli(count)=i;
        values(count)=1.0;
        count = count+1;
end

%Nodes
%  a-b-c-d-e
% nodec is center node

% Inner nodes
for i=3:sizex-2
        % Node indexes for computing finite differences
        % current node
        nodec = i;

        % surrounding nodes
        nodea = nodec-2;
        nodeb = nodec-1;
        noded = nodec+1;
        nodee = nodec+2;
        
        % Compute constants and derivatives in D
        A = D(nodec);
        B = (pm-pfill)*g;  
                      
        %Nodes
      %  a-b-c-d-e
    %k=  1-2-3-4-5
        
    %coefficients for 4th derivative, 2nd order accruacy
    % 1 -4 6 -4 1
               
     % indexes and values of design matrix
        for k=1:5
            rowi(count)=nodec;
            if k == 1
                coli(count)=nodea;
                values(count)=A/(delta^4.0);
            elseif k == 2
                coli(count)=nodeb;
                values(count)=-4*A/(delta^4.0);
            elseif k == 3
                coli(count)=nodec;
                values(count)=6*A/(delta^4.0)+B;
            elseif k == 4
                coli(count)=noded;
                values(count)=-4*A/(delta^4.0);
            elseif k == 5
                coli(count)=nodee;
                values(count)=A/(delta^4.0);
            
            end
            count = count+1;
        end
    end
    

% --------------------------------------------------
% Design Matrix (DM)
%size of sparse matrix is sizex by sizex
DM = sparse(rowi,coli,values,sizex,sizex);

% --------------------------------------------------
% Direct inversion to calculate deflections (w)
w = DM\tload;

% -------------------------------------------------
% Uncomment to write deflections out to txt file
% fid = fopen('deflection.txt','wt');
% for i=1:sizex
%         fprintf(fid,'%f\n',w(i));
% end        
% fclose(fid);



% Plot profiles 
figure(1)                % opening (and clearing) a new figure    
    subplot(4,1,[1 3]) 
    plot(x/1000,w,'-b',x/1000,-hload,'-r');  
    set(gca,'ydir','reverse')       % reversing the y-axis, so that the
                                    % y-axis increases downward
    %xlabel('Distance, km')    % adding in an x-axis label
    ylabel('Depth (m)')             % adding in an y-axis label
    %set(gca,'xaxislocation','top')  % putting the x-axis label on the top
                                    % of the plot
    box on;grid on; % grid minor;   % turning a plot border (box) on, and
                                    % overlying a simple grid
  %xlim([0 500])                  % setting the x-axis limits (deg C)
   % ylim([0 160])                   % setting the y-axis limits (km)
    tt1=title('Lithospheric Flexure Profile'); % adding a title to the plot
    set(tt1,'fontsize',12)           % making the title fontsize larger
    legend('Flexure','Load','Location','southeast') % add a legend inside plot   
    
    subplot(4,1,4)
    plot(x/1000,te/1000,'-k')
    xlabel('Distance, km')    % adding in an x-axis label
    ylabel('Te (km)')
    ylim([min(te)/1000-5 max(te)/1000+5])
    %set(gca,'yaxislocation','right')
    box on;grid on;
  
  
    