function out = lcf_strain_life_notch(mat_in,Kt,Sa,varargin)
%
%   lcf_strain_life_notch.m   ver 3.0   by Tom Irvine
%
%   Low cycle fatigue, strain-life method with Neuber notch analysis.
%
%   Worked example from:
%   https://blog.vibrationdata.com/2026/07/08/low-cycle-fatigue-strain-method/
%
%   ver 3.0 changes:
%     - Callable function with input arguments and an output struct
%     - Interactive prompts retained as the fallback for missing arguments
%
%----------------------------------------------------------------------
%   USAGE
%
%     out = lcf_strain_life_notch(mat,Kt,Sa)
%     out = lcf_strain_life_notch(mat,Kt,Sa,'Name',value,...)
%     out = lcf_strain_life_notch                 % fully interactive
%
%   INPUTS
%
%     mat   material specification, one of:
%             []                empty, prompts from the library
%             numeric scalar    index into the built-in library
%             char string       name or partial name match, e.g. '4340'
%             struct            user-supplied properties with fields
%                                 name, E, sfp, b, efp, c
%                               and optional fields
%                                 sy, np, Kp
%
%     Kt    elastic stress concentration factor, dimensionless
%     Sa    nominal stress amplitude, MPa
%
%   OPTIONAL NAME-VALUE PAIRS
%
%     'sm'       local mean stress, MPa.  Default 0.
%     'plots'    true or false, generate the three figures.  Default true.
%     'verbose'  true or false, print the report.  Default true.
%     'figstart' first figure number.  Default 1.
%
%   OUTPUT
%
%     out   struct with fields
%             mat        material properties actually used
%             Kt, Sa, sm inputs echoed
%             C          Neuber constant, MPa
%             sa, ea     local stress amplitude (MPa), strain amplitude
%             ea_el      elastic strain amplitude
%             ea_pl      plastic strain amplitude
%             dsig, deps stabilized loop ranges
%             twoNf, Nf  reversals and cycles to failure
%             term_el    elastic term of the strain-life equation
%             term_pl    plastic term
%             twoNt      transition life, reversals
%             twoNf_morrow, twoNf_swt   mean stress corrected lives
%             warn       cell array of warning messages, empty if none
%
%   EXAMPLES
%
%     out = lcf_strain_life_notch(1,3.0,300);
%     out = lcf_strain_life_notch('7075',2.5,180,'plots',false);
%     out = lcf_strain_life_notch([],3.0,300,'sm',50);
%
%     m.name='Custom'; m.E=2e5; m.sfp=1900; m.b=-0.09; m.efp=0.6; m.c=-0.6;
%     out = lcf_strain_life_notch(m,3.0,300);
%
%     Batch sweep over Kt:
%       Ktv = 1.5:0.25:4;
%       for i=1:length(Ktv)
%           r(i) = lcf_strain_life_notch(1,Ktv(i),300,'plots',false,...
%                                        'verbose',false);
%       end
%       loglog(Ktv,[r.Nf],'-o'); grid on;
%       xlabel('Kt'); ylabel('Cycles to Failure');
%
%   Units:  MPa, dimensionless strain
%
%----------------------------------------------------------------------

   if(nargin==0)
       fprintf('\n  lcf_strain_life_notch.m  ver 3.0  by Tom Irvine \n');
   end

%----------------------------------------------------------------------
%  Parse optional arguments
%----------------------------------------------------------------------

   opt = parse_options(varargin);

%----------------------------------------------------------------------
%  Resolve the material
%----------------------------------------------------------------------

   if(nargin < 1)
       mat_in = [];
   end

   mat = resolve_material(mat_in);
   mat = check_material(mat);

%  Cyclic deformation constants from the consistency relations
%     np = b/c          Kp = sfp/(efp^np)
%
%  A supplied np and Kp take precedence.  Otherwise derive them.

   if( ~isfield(mat,'np') || isempty(mat.np) )
       mat.np = mat.b/mat.c;
       mat.Kp = mat.sfp/(mat.efp^mat.np);
       mat.cyc_source = 'compatibility relations';
   else
       if( ~isfield(mat,'Kp') || isempty(mat.Kp) )
           mat.Kp = mat.sfp/(mat.efp^mat.np);
       end
       mat.cyc_source = 'supplied values';
   end

%----------------------------------------------------------------------
%  Resolve loading and geometry
%----------------------------------------------------------------------

   if(nargin < 2 || isempty(Kt))
       Kt = input_scalar('  Enter elastic stress concentration factor Kt',...
                         3.0, 1.0, 20.0);
   else
       validate_scalar(Kt,'Kt',1.0,20.0);
   end

   if(nargin < 3 || isempty(Sa))
       Sa = input_scalar('  Enter nominal stress amplitude Sa (MPa), fully reversed',...
                         300, 0, 1e5);
   else
       validate_scalar(Sa,'Sa',0,1e5);
   end

   sm = opt.sm;

   warn = {};

%  Sanity check.  The Neuber solution assumes the nominal section
%  remains essentially elastic.

   if( isfield(mat,'sy') && ~isempty(mat.sy) && Sa > mat.sy )
       warn{end+1} = sprintf(['nominal Sa = %.0f MPa exceeds the yield strength ',...
                              'of %.0f MPa; gross section yielding invalidates ',...
                              'the Neuber notch assumption'],Sa,mat.sy);
   end

   if( sm >= mat.sfp )
       error('lcf:meanStress',...
             'Local mean stress sm = %g MPa must be less than sf'' = %g MPa.',...
             sm,mat.sfp);
   end

%----------------------------------------------------------------------
%  Step 1.  Neuber constant
%----------------------------------------------------------------------

   C = (Kt*Sa)^2/mat.E;

%----------------------------------------------------------------------
%  Step 2.  Solve Neuber hyperbola with the cyclic stress-strain curve
%
%     sa*[ sa/E + (sa/Kp)^(1/np) ] = C
%
%  Solved by bisection so that no toolbox or fzero dependency is
%  required.  The residual is monotonic increasing in sa, so bisection
%  is robust.
%----------------------------------------------------------------------

   fneuber = @(s) s.*( s/mat.E + (s/mat.Kp).^(1/mat.np) ) - C;

   slo = 1e-6;
   shi = Kt*Sa;         % elastic estimate bounds the solution from above

   if(Sa == 0)
       sa = 0;  ea = 0;  ea_el = 0;  ea_pl = 0;
   else
       sa    = bisect(fneuber,slo,shi,1e-10,200);
       ea    = C/sa;              % total local strain amplitude
       ea_el = sa/mat.E;          % elastic part
       ea_pl = ea - ea_el;        % plastic part
   end

%----------------------------------------------------------------------
%  Step 3.  Stabilized hysteresis loop, Massing behavior
%----------------------------------------------------------------------

   dsig = 2*sa;
   deps = 2*ea;

%----------------------------------------------------------------------
%  Step 4.  Invert the strain-life equation
%
%     ea = (sfp/E)*(2Nf)^b + efp*(2Nf)^c
%
%  Monotonic decreasing in 2Nf, so bisection again.
%----------------------------------------------------------------------

   flife = @(x) (mat.sfp/mat.E)*x.^mat.b + mat.efp*x.^mat.c - ea;

   if( ea > 0 && flife(1) < 0 )

       warn{end+1} = ['local strain amplitude exceeds the strain-life curve ',...
                      'at 2Nf = 1; failure on the first reversal'];

       twoNf   = 1;
       Nf      = 0.5;
       term_el = (mat.sfp/mat.E);
       term_pl = mat.efp;

   elseif( ea == 0 )

       twoNf = Inf;  Nf = Inf;  term_el = 0;  term_pl = 0;

   else

       twoNf   = bisect(flife,1,1e12,1e-10,300);
       Nf      = twoNf/2;
       term_el = (mat.sfp/mat.E)*twoNf^mat.b;
       term_pl = mat.efp*twoNf^mat.c;

   end

%----------------------------------------------------------------------
%  Transition life
%----------------------------------------------------------------------

   twoNt = ( mat.efp*mat.E/mat.sfp )^( 1/(mat.b-mat.c) );

%----------------------------------------------------------------------
%  Mean stress corrections.  For sm = 0 these reduce to the
%  uncorrected solution.
%----------------------------------------------------------------------

   smax = sm + sa;

   if( ea > 0 )

       twoNf_morrow = bisect(@(x) ((mat.sfp-sm)/mat.E)*x.^mat.b ...
                                  + mat.efp*x.^mat.c - ea, 1,1e12,1e-10,300);

       if( smax > 0 )
           twoNf_swt = bisect(@(x) (mat.sfp^2/mat.E)*x.^(2*mat.b) ...
                                   + mat.sfp*mat.efp*x.^(mat.b+mat.c) ...
                                   - smax*ea, 1,1e12,1e-10,300);
       else
%          Compressive smax, SWT predicts no damage
           twoNf_swt = Inf;
       end

   else
       twoNf_morrow = Inf;
       twoNf_swt    = Inf;
   end

%----------------------------------------------------------------------
%  Assemble the output struct
%----------------------------------------------------------------------

   out.mat          = mat;
   out.Kt           = Kt;
   out.Sa           = Sa;
   out.sm           = sm;
   out.C            = C;
   out.sa           = sa;
   out.ea           = ea;
   out.ea_el        = ea_el;
   out.ea_pl        = ea_pl;
   out.dsig         = dsig;
   out.deps         = deps;
   out.twoNf        = twoNf;
   out.Nf           = Nf;
   out.term_el      = term_el;
   out.term_pl      = term_pl;
   out.twoNt        = twoNt;
   out.twoNf_morrow = twoNf_morrow;
   out.twoNf_swt    = twoNf_swt;
   out.warn         = warn;

%----------------------------------------------------------------------
%  Report and plots
%----------------------------------------------------------------------

   if(opt.verbose)
       print_report(out);
   end

   if(opt.plots && ea > 0)
       make_plots(out,opt.figstart);
   end

end
%======================================================================
function print_report(r)
%
%  Print the formatted analysis report.
%
   mat = r.mat;

   fprintf('\n  Material:  %s \n',mat.name);
   fprintf('    E    = %8.0f MPa \n',mat.E);
   fprintf('    sf''  = %8.0f MPa \n',mat.sfp);
   fprintf('    b    = %8.3f \n',mat.b);
   fprintf('    ef''  = %8.3f \n',mat.efp);
   fprintf('    c    = %8.3f \n',mat.c);
   fprintf('    n''   = %8.4f   (%s) \n',mat.np,mat.cyc_source);
   fprintf('    K''   = %8.0f MPa \n',mat.Kp);

   fprintf('\n  Loading: \n');
   fprintf('    Kt   = %8.2f \n',r.Kt);
   fprintf('    Sa   = %8.1f MPa \n',r.Sa);
   fprintf('    sm   = %8.1f MPa   (local mean stress) \n',r.sm);

   for i=1:length(r.warn)
       fprintf('\n  WARNING: %s. \n',r.warn{i});
   end

   fprintf('\n  Step 1.  Neuber constant \n');
   fprintf('    sa*ea = (Kt*Sa)^2/E = %8.4f MPa \n',r.C);

   fprintf('\n  Step 2.  Local notch stress and strain \n');
   fprintf('    sigma_a       = %8.1f MPa \n',r.sa);
   fprintf('    epsilon_a     = %11.4e \n',r.ea);
   fprintf('      elastic     = %11.4e \n',r.ea_el);
   fprintf('      plastic     = %11.4e \n',r.ea_pl);
   if(r.ea > 0)
       fprintf('      plastic frac= %8.1f %% \n',100*r.ea_pl/r.ea);
   end

   fprintf('\n    Elastic estimate vs Neuber solution: \n');
   fprintf('      stress:  %7.1f MPa  ->  %7.1f MPa   (elastic overestimates) \n',...
                  r.Kt*r.Sa,r.sa);
   fprintf('      strain: %10.4e  -> %10.4e   (elastic underestimates) \n',...
                  r.Kt*r.Sa/mat.E,r.ea);

   fprintf('\n  Step 3.  Stabilized hysteresis loop (Massing) \n');
   fprintf('    delta_sigma   = %8.0f MPa \n',r.dsig);
   fprintf('    delta_epsilon = %11.4e \n',r.deps);

   fprintf('\n  Step 4.  Strain-life solution \n');
   fprintf('    2Nf           = %10.0f reversals \n',r.twoNf);
   fprintf('    Nf            = %10.0f cycles \n',r.Nf);
   fprintf('      elastic term= %11.4e \n',r.term_el);
   fprintf('      plastic term= %11.4e \n',r.term_pl);

   fprintf('\n  Transition life \n');
   fprintf('    2Nt           = %10.0f reversals \n',r.twoNt);

   if(r.twoNf > r.twoNt)
       fprintf('    Predicted life is ABOVE the transition life. \n');
       fprintf('    Mixed elastic-plastic zone; plastic strain still \n');
       fprintf('    carries a meaningful share of the damage. \n');
   else
       fprintf('    Predicted life is BELOW the transition life. \n');
       fprintf('    Plastic strain dominates; ductility governs. \n');
   end

   fprintf('\n  Mean stress corrections, sigma_m = %6.1f MPa \n',r.sm);
   fprintf('    Morrow  2Nf   = %10.0f reversals \n',r.twoNf_morrow);
   fprintf('    SWT     2Nf   = %10.0f reversals \n',r.twoNf_swt);

   if(r.sm == 0)
       fprintf('    Both reduce to the uncorrected solution for sm = 0. \n');
   end

   fprintf('\n  Done. \n\n');

end
%======================================================================
function make_plots(r,fs)
%
%  Generate the three figures.
%
   mat = r.mat;
   Kt  = r.Kt;
   Sa  = r.Sa;
   sa  = r.sa;
   ea  = r.ea;

%  --- Plot 1.  Neuber hyperbola and the cyclic stress-strain curve ---

   svec = linspace(1,1.05*Kt*Sa,600);
   ecyc = svec/mat.E + (svec/mat.Kp).^(1/mat.np);   % cyclic curve
   ehyp = r.C./svec;                                % Neuber hyperbola

   figure(fs);
   plot(1000*ecyc,svec,'k-','LineWidth',1.5); hold on;
   plot(1000*ehyp,svec,'r--','LineWidth',1.5);
   plot(1000*ea,sa,'bo','MarkerSize',9,'MarkerFaceColor','b');
   grid on;
   xlabel('Strain Amplitude (millistrain)');
   ylabel('Stress Amplitude (MPa)');
   title(['Neuber Hyperbola and Cyclic Stress-Strain Curve, ',mat.name]);
   legend('Cyclic stress-strain curve','Neuber hyperbola',...
          'Local notch solution','Location','SouthEast');
   xlim([0 1.6*1000*ea]);
   ylim([0 1.05*Kt*Sa]);

%  --- Plot 2.  Stabilized hysteresis loop, Massing ---

   ds = linspace(0,r.dsig,400);
   de = ds/mat.E + 2*(ds/(2*mat.Kp)).^(1/mat.np);

%  Upper branch, from lower tip upward.  Lower branch by point symmetry.

   e_up =  -ea + de;      s_up = -sa + ds;
   e_dn =   ea - de;      s_dn =  sa - ds;

   figure(fs+1);
   plot(1000*e_up,s_up,'b-','LineWidth',1.5); hold on;
   plot(1000*e_dn,s_dn,'b-','LineWidth',1.5);
   plot(1000*[-ea ea],[-sa sa],'ro','MarkerSize',8,'MarkerFaceColor','r');
   grid on;
   xlabel('Local Strain (millistrain)');
   ylabel('Local Stress (MPa)');
   title(['Stabilized Hysteresis Loop, Massing Behavior, ',mat.name]);
   legend('Loop branches','Loop tips','Location','SouthEast');

%  --- Plot 3.  Strain-life curve with the solution marked ---

   x  = logspace(0,8,600);
   te = (mat.sfp/mat.E)*x.^mat.b;
   tp = mat.efp*x.^mat.c;

   figure(fs+2);
   loglog(x,te+tp,'k-','LineWidth',1.8); hold on;
   loglog(x,te,'b--','LineWidth',1.2);
   loglog(x,tp,'r--','LineWidth',1.2);
   loglog(r.twoNf,ea,'ko','MarkerSize',9,'MarkerFaceColor','g');
   loglog([r.twoNt r.twoNt],[1e-5 1],'m:','LineWidth',1.5);
   grid on;
   xlabel('Reversals to Failure, 2Nf');
   ylabel('Strain Amplitude');
   title(['Strain-Life Curve, ',mat.name]);
   legend('Total strain','Elastic term','Plastic term',...
          'Solution','Transition life','Location','SouthWest');
   ylim([1e-4 1e-1]);

end
%======================================================================
function opt = parse_options(args)
%
%  Parse name-value pairs with defaults.
%
   opt.sm       = 0;
   opt.plots    = true;
   opt.verbose  = true;
   opt.figstart = 1;

   if( mod(length(args),2) ~= 0 )
       error('lcf:options','Optional arguments must be Name,value pairs.');
   end

   for i=1:2:length(args)

       nm = args{i};
       vl = args{i+1};

       if(~ischar(nm))
           error('lcf:options','Option name must be a string.');
       end

       switch lower(nm)
           case 'sm'
               validate_scalar(vl,'sm',-1e5,1e5);
               opt.sm = vl;
           case 'plots'
               opt.plots = logical(vl);
           case 'verbose'
               opt.verbose = logical(vl);
           case 'figstart'
               validate_scalar(vl,'figstart',1,1e4);
               opt.figstart = round(vl);
           otherwise
               error('lcf:options','Unrecognized option: %s',nm);
       end

   end

end
%======================================================================
function mat = resolve_material(m)
%
%  Resolve the material argument into a property struct.
%
   lib = material_library;

   if( isempty(m) )
       mat = select_material(lib);

   elseif( isstruct(m) )
       mat = m;

   elseif( isnumeric(m) && isscalar(m) )
       k = round(m);
       if( k < 1 || k > length(lib) )
           error('lcf:matIndex',...
                 'Material index %d out of range, library has %d entries.',...
                 k,length(lib));
       end
       mat = lib(k);

   elseif( ischar(m) )
       names = {lib.name};
       hit   = find( ~cellfun(@isempty, strfind(lower(names),lower(m))) );

       if( isempty(hit) )
           error('lcf:matName','No library material matches "%s".',m);
       end

       if( length(hit) > 1 )
           fprintf('\n  "%s" matches more than one library entry: \n',m);
           for i=1:length(hit)
               fprintf('    %2d  %s \n',hit(i),lib(hit(i)).name);
           end
           error('lcf:matName','Ambiguous material name, be more specific.');
       end

       mat = lib(hit);

   else
       error('lcf:matType',...
             'Material must be empty, an index, a name string, or a struct.');
   end

end
%======================================================================
function mat = check_material(mat)
%
%  Verify that a supplied material struct has the required fields
%  with physically sensible values.
%
   req = {'E','sfp','b','efp','c'};

   for i=1:length(req)
       if( ~isfield(mat,req{i}) || isempty(mat.(req{i})) )
           error('lcf:matField','Material is missing required field "%s".',req{i});
       end
       v = mat.(req{i});
       if( ~isnumeric(v) || ~isscalar(v) || ~isreal(v) || ~isfinite(v) )
           error('lcf:matField','Material field "%s" must be a finite real scalar.',req{i});
       end
   end

   if( ~isfield(mat,'name') || isempty(mat.name) )
       mat.name = 'User-defined material';
   end

   if( mat.E   <= 0 ), error('lcf:matField','E must be positive.');   end
   if( mat.sfp <= 0 ), error('lcf:matField','sf'' must be positive.'); end
   if( mat.efp <= 0 ), error('lcf:matField','ef'' must be positive.'); end
   if( mat.b   >= 0 ), error('lcf:matField','b must be negative.');   end
   if( mat.c   >= 0 ), error('lcf:matField','c must be negative.');   end

%  The transition life exponent 1/(b-c) requires b ~= c

   if( abs(mat.b - mat.c) < 1e-9 )
       error('lcf:matField','b and c must differ; transition life is undefined.');
   end

   if( ~isfield(mat,'sy') ), mat.sy = []; end
   if( ~isfield(mat,'np') ), mat.np = []; end
   if( ~isfield(mat,'Kp') ), mat.Kp = []; end

end
%======================================================================
function validate_scalar(v,nm,vmin,vmax)
%
%  Verify an argument is a finite real scalar within range.
%
   if( ~isnumeric(v) || ~isscalar(v) || ~isreal(v) || ~isfinite(v) )
       error('lcf:arg','%s must be a finite real scalar.',nm);
   end

   if( v < vmin || v > vmax )
       error('lcf:arg','%s = %g is outside the range [%g, %g].',nm,v,vmin,vmax);
   end

end
%======================================================================
function mat = select_material(lib)
%
%  Present the material library and return the selected properties.
%
   n = length(lib);

   fprintf('\n----------------------------------------------------------\n');
   fprintf('  Material Library \n');
   fprintf('----------------------------------------------------------\n');
   fprintf('   #                                E(MPa)   sf''(MPa)      b     ef''      c \n');
   fprintf('  ------------------------------------------------------------------------- \n');

   for i=1:n
       fprintf('  %2d  %-28s %7.0f %8.0f %7.3f %7.3f %7.3f \n',...
               i,lib(i).name,lib(i).E,lib(i).sfp,lib(i).b,lib(i).efp,lib(i).c);
   end

   fprintf('  %2d  %-28s \n',n+1,'Enter properties manually');
   fprintf('  ------------------------------------------------------------------------- \n');

   k = round( input_scalar('\n  Select material number',1,1,n+1) );

   if(k <= n)
       mat = lib(k);
   else
       mat = manual_material;
   end

end
%======================================================================
function mat = manual_material
%
%  Prompt for a user-defined material.
%
   fprintf('\n  Enter material properties. \n');

   nm = input('  Material name: ','s');
   if(isempty(nm))
       nm = 'User-defined material';
   end

   mat.name = nm;
   mat.E    = input_scalar('  Elastic modulus E (MPa)',           200000, 1,   1e7);
   mat.sfp  = input_scalar('  Fatigue strength coefficient sf'' (MPa)',1900, 1,   1e5);
   mat.b    = input_scalar('  Fatigue strength exponent b',        -0.09, -1,   -0.01);
   mat.efp  = input_scalar('  Fatigue ductility coefficient ef''',    0.6, 1e-4, 10);
   mat.c    = input_scalar('  Fatigue ductility exponent c',         -0.6, -2,   -0.1);
   mat.sy   = [];
   mat.np   = [];
   mat.Kp   = [];

end
%======================================================================
function lib = material_library
%
%  Built-in strain-life material library.
%
%  Fields:
%    name   description including heat treatment or temper
%    E      elastic modulus, MPa
%    sy     approximate monotonic yield strength, MPa, for a warning check
%    sfp    fatigue strength coefficient sf', MPa
%    b      fatigue strength exponent
%    efp    fatigue ductility coefficient ef'
%    c      fatigue ductility exponent
%    np     cyclic strain hardening exponent n', empty to derive as b/c
%    Kp     cyclic strength coefficient K', MPa, empty to derive
%
%  Values are representative textbook figures for demonstration.
%  Verify against the governing material specification or test data
%  before use in a design analysis.
%
   k = 0;

%  --- Steels ---

   k = k+1;
   lib(k) = mk('SAE 4340 steel, Q&T',      200000, 1172, 1900, -0.09,  0.60, -0.60, [],   []);

   k = k+1;
   lib(k) = mk('SAE 4340 steel, 350 HB',   200000, 1240, 1758, -0.076, 0.45, -0.62, 0.14, 1940);

   k = k+1;
   lib(k) = mk('SAE 1045 steel, 225 HB',   200000,  634, 1225, -0.095, 1.00, -0.66, 0.18, 1250);

   k = k+1;
   lib(k) = mk('SAE 1045 steel, 500 HB',   200000, 1720, 2725, -0.081, 0.07, -0.60, 0.12, 3330);

   k = k+1;
   lib(k) = mk('SAE 1020 steel, hot rolled',203000, 262,  895, -0.12,  0.41, -0.51, 0.18,  941);

   k = k+1;
   lib(k) = mk('ASTM A36 structural steel', 200000, 250,  717, -0.104, 0.27, -0.48, 0.22,  976);

   k = k+1;
   lib(k) = mk('AISI 304 stainless, ann.',  193000, 240, 1000, -0.114, 0.17, -0.402,0.28, 1660);

   k = k+1;
   lib(k) = mk('AISI 316 stainless, ann.',  193000, 250, 1000, -0.11,  0.32, -0.47, 0.27, 1600);

   k = k+1;
   lib(k) = mk('Maraging 300 steel',        186000, 1930, 2400, -0.09, 0.30, -0.70, [],   []);

%  --- Aluminum alloys ---

   k = k+1;
   lib(k) = mk('2024-T3 aluminum',           73100,  345, 1100, -0.124, 0.22, -0.59, 0.20,  807);

   k = k+1;
   lib(k) = mk('2024-T351 aluminum',         73100,  380, 1100, -0.124, 0.22, -0.59, 0.20,  807);

   k = k+1;
   lib(k) = mk('7075-T6 aluminum',           71000,  469, 1466, -0.143, 0.262,-0.619,0.146, 913);

   k = k+1;
   lib(k) = mk('6061-T6 aluminum',           68900,  276,  545, -0.093, 0.35, -0.60, 0.14,  462);

   k = k+1;
   lib(k) = mk('5083-H116 aluminum',         71000,  228,  540, -0.09,  0.40, -0.63, [],   []);

%  --- Titanium ---

   k = k+1;
   lib(k) = mk('Ti-6Al-4V, annealed',       117000,  930, 2030, -0.104, 0.841,-0.688,0.10, 1772);

%  --- Nickel alloys ---

   k = k+1;
   lib(k) = mk('Inconel 718, aged',         200000, 1100, 1640, -0.06,  2.67, -0.82, 0.08, 1910);

   k = k+1;
   lib(k) = mk('Haynes 188, annealed',      206000,  464, 1150, -0.08,  0.35, -0.60, [],   []);

%  --- Copper alloys ---

   k = k+1;
   lib(k) = mk('OFHC copper, annealed',     115000,   70,  400, -0.11,  0.55, -0.62, [],   []);

end
%======================================================================
function s = mk(name,E,sy,sfp,b,efp,c,np,Kp)
%
%  Assemble one material library entry.
%
   s.name = name;
   s.E    = E;
   s.sy   = sy;
   s.sfp  = sfp;
   s.b    = b;
   s.efp  = efp;
   s.c    = c;
   s.np   = np;
   s.Kp   = Kp;
end
%======================================================================
function v = input_scalar(prompt,default,vmin,vmax)
%
%  Prompt for a scalar with a default value and range checking.
%  A blank entry accepts the default.
%
   while(1)

       str = input(sprintf('%s [%g]: ',prompt,default),'s');

       if(isempty(str))
           v = default;
           return;
       end

       v = str2double(str);

       if( isnan(v) || ~isreal(v) )
           fprintf('    Invalid entry.  Enter a number. \n');
           continue;
       end

       if( v < vmin || v > vmax )
           fprintf('    Out of range.  Enter a value between %g and %g. \n',vmin,vmax);
           continue;
       end

       return;

   end
end
%======================================================================
function x = bisect(f,a,b,tol,itmax)
%
%  Simple bisection root finder.  Requires a sign change on [a,b].
%
   fa = f(a);
   fb = f(b);

   if(fa*fb > 0)
       error('bisect: no sign change on the given interval');
   end

   for i=1:itmax
       x  = 0.5*(a+b);
       fx = f(x);
       if( abs(fx) < tol || 0.5*(b-a) < tol*max(1,abs(x)) )
           return;
       end
       if(fa*fx < 0)
           b = x;  fb = fx;
       else
           a = x;  fa = fx;
       end
   end
end