

%lab 5 - lecture notes chapter 3 - plotting the graph of an inverse function
graphics_toolkit("gnuplot");
disp("Sigma function");

function v=Sigma(s,M=0.5)
%Sigma function M\in[0,1)
%s\in[0,1]

v=(2/(M-1))*(1-(1-s).^(M-1)) -s.*(1-s).^(M-1);

end  

%first plot of the graph
s=linspace(0,0.99,300);
y=Sigma(s);
plot(s,y,";Sigma;");
pause(1);
plot(y,s,"+;Sigma inverse;");
pause(1);

%inverse of sigma

function [s,fc]=sigmainv(z,M=0.5)
%computes the inverse of Sigma for z

 if(z==0)
  s=0;fc=0;
 else
   [s,FS,info,out]=fzero(@(s) Sigma(s,M)-z,[0,1-eps/2]);
   if(info!=1)
    warning("Something wrong with fzero()");
   end
   fc=out.funcCount; 
 end
end 

M=0.5;
y=linspace(0,60,300);
s=fc=zeros(size(y));
for k=1:length(y),
  [s(k),fc(k)]=sigmainv(y(k),M);
end  
plot(y,s,"+;solve S(s)=z;");
pause(1);
allfc=sum(fc)

function v=G(s,z,M=0.5)
  p=2/(1-M);
  v=1-s -((p-s)./(p+z)).^(p/2);
end  
s=rand;z=Sigma(s);G(s,z)
 
function [s,fc]=sigmainv1(z,M=0.5)
%computes the inverse of Sigma for z

 if(z==0)
  s=0;fc=0;
 else

  [s,FS,info,out]=fzero(@(s) G(s,z,M),[0.1-eps/2]);
   if(info!=1)
    warning("Something wrong with fzero()");
   end
   fc=out.funcCount; 
 end
end 

M=0.5;
y=linspace(0,60,300);
s1=fc1=zeros(size(y));
for k=1:length(y),
  [s1(k),fc(k)]=sigmainv1(y(k),M);
end  
plot(y,s1,"+;Gz;");
pause(1);
allfc1=sum(fc1)

%we create a very very dense mesh getting denser close to onenormest
s=1-logspace(-13,0,10000);
y=Sigma(s);
good=(y<=60);
y=y(good);
s=s(good);
plot(y,s,"+;sigma inv;");



