修正剑桥模型程序

%% MATLAB Code for Cam-Clay Model 
clear all;
%% License
out1=fprintf（‘\n \t\t   MATLAB CODE FOR SIMULATION OF MODIFIED CAMCLAY\n‘）;
out2=fprintf（‘\t Copyright （C） 2011 Krishna Kumar University of Cambridge\n‘）;
out11=fprintf（‘\n\t\t\t The program is distributed under GNU GPL v 2.0 ‘）;
out12=fprintf（‘\n\t\t view license agreement at http://www.gnu.org/licenses/\n‘）;
%{   
    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation either version 3 of the License or
    （at your option） any later version.

    This program is distributed in the hope that it will be useful
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program.  If not see ‘）;
%}
%% Input Parameters
    out3=fprintf（‘\nINPUT PARAMETERS FOR MODIFIED CAMCLAY\n\n‘）;
    cp=input（‘Enter the inital Consolidation pressure （kPa） （eg. 150 kPa）  = ‘）;%cp=150;   
    p0=input（‘Enter the initial Confining pressure （kPa）    （eg. 150 kPa）  = ‘）;%p0=150;
    M=input（‘Enter the value of Critical Friction Anlge M   （eg. 0.95） = ‘）;%M=0.95;
    l=input（‘Enter the value of Lamda                       （eg. 0.2） = ‘）;%l=0.2;
    k=input（‘Enter the value of Kappa                       （eg. 0.04） = ‘）;%k=0.04;    
    N=input（‘Enter the value of N                           （eg. 2.5） = ‘）;%N=2.5;
    v=input（‘Enter the value of poissons ratio              （eg. 0.15） = ‘）;%v=0.15;
%% Computation of Other Parameters （Ve0 and OCR）
    pc=cp;V=N-（l*log（pc））;e0=V-1;OCR=cp/p0;%Initalizing confining pressure
%% Strain Increament and Strain Matrix Definition
    out4=fprintf（‘\nSTRAIN INCREAMENT AND ITERATION\n\n‘）;
    iteration=input（‘Enter number of iterations to perform  （eg. 5000） = ‘）;
    if （iteration <=2000 || iteration == ‘ ‘）
    iter=2000;
    out5=fprintf（‘\nThe iterations entered is too low using defaults %f\n‘iter）;
    else iter=iteration;
    end
    strsteps=input（‘Enter the strain increament （in decimal） （eg. 0.01） = ‘）;
    if （strsteps > 0.01|| strsteps == ‘ ‘ || strsteps <= 0.）
    ide=0.01;
    out6=fprintf（‘\nThe strain step entered is too low using defaults %f\n‘ide）;
    else ide=strsteps;
    end
    de=ide/100;
    es=0:ide:（iter-1）*ide; %strain
    dstrain=[de;-de/2;-de/2;0;0;0];%strain increament
%% Block Memory allocation
    De=zeros（16）;
    dfds=zeros（61）;
    dfdep=zeros（61）;
    u=zeros（1iter）;
    p=zeros（1iter）;
    q=zeros（1iter）;
%% Yield Surface and Conditions
    p1=（0:pc）;% CSL in p-q space
    q1 = M*p1;
    qy=（M^2*（pc*p1-p1.^2））.^0.5;%Plot the initial yield locus
%% Initialize   
     a=1;
     S=[p0;p0;p0;0;0;0];