clear all

  e = 0.001;
  et = 0.1;             
  E_nom = 100;                 % E_nom = 100e+9 (Pa)
  c_nom = 0.1;                 % c_nom = 0.1 (Pa s)
  s = E_nom*e + c_nom*et;

  per = 0.1;
  npts = 200;
  Ept1 = E_nom - per*E_nom;
  Ept2 = E_nom + per*E_nom;
  cpt1 = c_nom - per*c_nom;
  cpt2 = c_nom + per*c_nom;
  Edt = (Ept2-Ept1)/npts;  
  cdt = (cpt2-cpt1)/npts;

  ypt1 = 0.095;
  ypt2 = .125;

  Evals = Ept1:Edt:Ept2;
  cvals = cpt1:cdt:cpt2; 
  Emean = e*Evals + et*c_nom; 
  cmean = et*cvals + e*E_nom;  

  var1 = ((Ept2-Ept1)^2)/12;
  var2 = ((cpt2-cpt1)^2)/12;
  denom = e^2*var1 + et^2*var2;

  figure(1)
  plot([Ept1 Ept2],[s s],'k','linewidth',2)
  axis([Ept1 Ept2 ypt1 ypt2])
  hold on
  plot([100 100],[ypt1 ypt2],'k','linewidth',2)
  plot(Evals,Emean,'-b','linewidth',3)
  hold off
  set(gca,'Fontsize',[22]);
  box on
  xlabel('E (GPa)')
  ylabel('s (GPa)')

  figure(2)
  plot([cpt1 cpt2],[s s],'k','linewidth',2)
  axis([cpt1 cpt2 ypt1 ypt2])
  hold on
  plot([.1 .1],[ypt1 ypt2],'k','linewidth',2)
  plot(cvals,cmean,'-b','linewidth',3)
  hold off
  set(gca,'Fontsize',[22]);
  box on
  xlabel('c (GPa\cdots)')
  ylabel('s (GPa)')