disp(" ")

%
% SECTION 1: Lines y = F(x).         
%

%
% Cubic line
%

function y = Cube(x)
a = 1; b = 0; c = -1; d = 0; 	% Coefficients       
y = a*x.^3 + b*x.^2 + c*x + d;
endfunction

%
% SECTION 2: Lines (x,y) = F(t). 
%

%
% Circle
%

function [x,y] = Circle(t)
a = 1; 				% Radius
x = a*sin(t);               
y = a*cos(t);               
endfunction           

%
% Archimedian Spiral      
%

function [x,y] = Arch(t)
R = 1; 				% Radius
a = 1;				% Angular speed
x = R*t*sin(a*t);
y = R*t*cos(a*t);
endfunction;


%
% SECTION 3: Lines (x,y,z) = F(t). 
%

%
% Elicoid
%

function [x,y,z] = Elicoid(t)
R = 1;                        	% Radius
a = 1;                        	% Angular Speed    
v = 1;                        	% Vertical Speed      
x = R*sin(a*t);               
y = R*cos(a*t);               
z = v*t;               
endfunction           

%
% Tornado
%

function [x,y,z] = Tornado(t)
x = t*sin(t);               
y = t*cos(t);               
z = t;               
endfunction           


%
% SECTION 4: Surfaces z = F(x,y). 
%

%
% Iperbolic Paraboloid
%

function z = IperParab(x,y)
a = 1;				% Axis a       
b = 1; 				% Axis b 
z = x.^2/a.^2 - y.^2/b.^2;               
endfunction           

%
% Parab
%

function z = Parab(x,y)
a = 1;				% Axis a       
b = 1; 				% Axis b 
z = x.^2/a.^2 + y.^2/b.^2;               
endfunction           

%
% Wave out of (a,b)       
%

function z = W(a,b,x,y)
z = sin(sqrt((x-a).^2 + (y-b).^2));
endfunction

%
% Wave1              
%

function z = Wave1(x,y)
z = W(0,0,x,y);           
endfunction

%
% Wave2              
%

function z = Wave2(x,y)
z = W(0,0) + W(5,5);
endfunction

%
% Wave3                       
%

function z = Wave3(x,y)
z = W(0,0) + W(5,5) + W(5,0);
endfunction


%
% SECTION 5: Surfaces [x,y,z] = F(u,v). 
%

%
% Plane parametric on x, y coodinates   
%
 
function [x,y,z] = Plane(u,v)
x = u;               
y = v;               
z = 0;               
endfunction           

%
% Sphere 
%
 
function [x,y,z] = Sphere(u,v)
R = 1;				% Radius      
x = R*cos(u)*cos(v);               
y = R*sin(u);               
z = R*cos(u)*sin(v);               
endfunction           

%
% Cylinder
%

function [x,y,z] = Cylinder(u,v)
a = 1;   			% Radius
x = a*sin(v);               
y = a*cos(v);               
z = u;               
endfunction           


%
% Torus 
%

function [x,y,z] = Torus(a,b)
R = 1; r = 1;       		% Radius
x = cos(a)*(R+r*cos(b));
y = sin(a)*(R+r*cos(b));
z = r*sin(b);           
endfunction

%
% SECTION 6. Here starts an interface to Octave plotting tools. 
%

disp(" GUIDED TOUR: OCTAVE PLOTTING TOOLS."),       
disp("Course: Computer Tools for Maths - S. Berardi - Torino 1999.")
disp(" ")

function [a,b] = InputSegment(var1,var2)
 a = 0; b = 0;
 while (a >= b) 
  a = input(strcat(var1, " = "));           
  b = input(strcat(var2, " = "));      
  if (a >= b)
   m1 = strcat("Error: ",var1," >= ",var2,".");
   m2 = strcat(" [",var1,",",var2,"] is not a proper segment.");  
   disp(strcat(m1,m2))
  endif
 endwhile 
endfunction

function n = InputPositive(var)
 n = 0;
 while (n < 1)
  n = input(strcat(var, " = "));
  if (n < 1)
   m1 = strcat("Error: ",var," < 1.");
   m2 = strcat(" It is not correct to divide a segment by ", var);  
  endif
 endwhile
endfunction

%
% Line y = Y(x)
%

function g = LineDraw(Y,x1,x2,n)
 s = (x2 - x1)/n;
  for i = 1:(n+1)
   g(i,1) = x1 + (i-1)*s;
   g(i,2) = feval(Y,g(i,1));
  endfor
 endfunction


function LinePlot(Y)

 [x1,x2] = InputSegment("x1","x2");
 n       = InputPositive("n");

 g = zeros(n+1,2);         
 
 M =  LineDraw(Y,x1,x2,n);  
 gplot(M);
 save Plotting.data Y x1 x2 n M;
 disp("Data temporary saved in Plotting.data.")
 
endfunction


%
% Line (x,y) = F(t)        
%


function g = LineParDraw(F,t1,t2,n)
 s = (t2 - t1)/n;
 for i = 1:(n+1)
  time   = t1 + (i-1)*s;
  [g(i,1),g(i,2)] = feval(F,time);
 endfor
endfunction

function LineParPlot(F)      

 [t1,t2] = InputSegment("t1","t2");
 n       = InputPositive("n");

 g = zeros(n+1,2);         
 
 M =  LineParDraw(F,t1,t2,n);  
 disp("Plot starting at time t1, ending at time t2, step for time:"), 
 disp((t2 - t1)/n), 
 gplot(M);
 save Plotting.data F t1 t2 n M;
 disp("Data temporary saved in Plotting.data.") 
endfunction


%
% Line 3D parametric
%


function g = Line3DDraw(F,t1,t2,n)
s = (t2 - t1)/n;
 for i = 1:(n+1)
  time   = t1 + (i-1)*s;
  [g(i,1),g(i,2),g(i,3)] = feval(F,time);
 endfor
endfunction

function Line3DPlot(F)       

 gset parametric

 [t1,t2] = InputSegment("t1","t2");
 n       = InputPositive("n");
 
 g = zeros(n+1,3);         
 
 M =  Line3DDraw(F,t1,t2,n);  
 disp("Plot starting at time t1, ending at time t2, step for time:"), 
 disp((t2 - t1)/n), 
 gsplot(M);
 save Plotting.data F t1 t2 n M ;
 disp("Data temporary saved in Plotting.data")
endfunction

 
% We define an (n+1)x(n+1) matrix storing all required values of f. 


function g = SurfaceDraw(Z,x1,y1,x2,y2,n)
s = (x2 - x1)/n;
t = (y2 - y1)/n;
 for i = 1:(n+1)
  for j = 1:(n+1)
   g(i,j) = feval(Z,x1 + (i-1)*s, y1 + (j-1)*t);
  endfor
 endfor
endfunction

function SurfacePlot(Z)

 [x1,x2] = InputSegment("x1","x2");
 [y1,y2] = InputSegment("y1","y2");
 n       = InputPositive("n");

 g = zeros(n+1,n+1);         

 M =  SurfaceDraw(Z,x1,y1,x2,y2,n);  
 save Plotting.data Z x1 x2 y1 y2 n M;
 disp("Data temporary saved in Plotting.data.") 

 axis();
 mesh((x1:(x2-x1)/n:x2),(y1:(y2-y1)/n:y2),M);           

endfunction
     

%
% Parametric Surface
%


% We define a k x 2 matrix storing all required values, with k = 2*(n+1)*(m+1). 

function g = SurfaceParDraw(F,u1,v1,u2,v2,n,m)

s = (u2 - u1)/n;
t = (v2 - v1)/m;          
cursor = 1;
i = 1;

 while (i<= (n+1))
  for j = 1:(m+1)
   U           = u1 + (i-1)*s;
   V           = v1 + (j-1)*t;
   [g(cursor,1),g(cursor,2),g(cursor,3)] = feval(F,U,V);
   ++cursor; 
  endfor
  ++i;
  if (i <= (n+1))
   for j = 1:(m+1)
    U           = u1 + (i-1)*s;
    V           = v2 - (j-1)*t;
    [g(cursor,1),g(cursor,2),g(cursor,3)] = feval(F,U,V);
    ++cursor;
   endfor
  endif;
  ++i;
 endwhile               

 j = 1;

 while (j <= (m+1))
  for i = 1:(n+1)
   U           = u2 - (i-1)*s;
   V           = v2 - (j-1)*t;
   [g(cursor,1),g(cursor,2),g(cursor,3)] = feval(F,U,V);
   ++cursor; 
  endfor
  ++j;           
  if (j <= (m+1))
   for i = 1:(n+1)
    U           = u1 + (i-1)*s;
    V           = v2 - (j-1)*t;
    [g(cursor,1),g(cursor,2),g(cursor,3)] = feval(F,U,V);
    ++cursor;
   endfor
  endif
  ++j;
 endwhile               
 
endfunction


function SurfaceParPlot(F)

 [u1,u2] = InputSegment("u1","u2");
 [v1,v2] = InputSegment("v1","v2");
 n       = InputPositive("n");
 m       = InputPositive("m");

 g = zeros(2*(n+1)*(m+1),3);         
 
 M =  SurfaceParDraw(F,u1,v1,u2,v2,n,m);  
 disp("Plot parametrized over [u1,u2] X [v1,v2], steps for u, v:"), 
 disp((u2 - u1)/n), 
 disp((v2 - v1)/m),
 gset parametric
 gsplot(M)
 save Plotting.data F u1 v1 u2 v2 n m M;
 disp("Data temporary saved in Plotting.data.")

endfunction


function SavePlot(k)

file     = input("File name (without the extension .ps)", "s");
fileps   = strcat("\"",file,".ps\"");
filedata = strcat(file,".data");

load -force Plotting.data;
if (k == 1) 
 save(filedata, "Y", "x1", "x2", "n");                              endif;
if (k == 2) 
 save(filedata, "F",  "t1", "t2", "n");                        endif;
if (k == 3) 
 save(filedata, "F",  "Z",  "t1", "t2", "n");                  endif;
if (k == 4) 
 save(filedata, "Z", "x1", "x2", "y1", "y2", "n");                  endif;
if (k == 5) 
 save(filedata, "F",  "u1", "u2", "v1", "v2", "n", "m"); endif;

gset terminal postscript;
gset("out",fileps); % gset(.,.) requires input format: "X", "\"Y\"";
replot;  
gset terminal x11;
gset out "STDOUT";

disp(strcat("Data saved in ", filedata, ". Picture saved in ",fileps, "."))     
 
endfunction


function ChangePlot(k)

load -force Plotting.data; 

if (k == 1) 
 [y1,y2] = InputSegment("y1", "y2");
 axis([x1,x2,y1,y2]);
 replot;
endif;

if (k == 2) 
 [x1,x2] = InputSegment("x1", "x2");
 [y1,y2] = InputSegment("y1", "y2");
 axis([x1,x2,y1,y2]);
 replot;
endif;

if ((k == 3) | (k == 5))
 [x1,x2] = InputSegment("x1", "x2");
 [y1,y2] = InputSegment("y1", "y2");
 [z1,z2] = InputSegment("z1", "z2");
 axis([x1,x2,y1,y2,z1,z2]);
 replot;
endif;

if (k == 4) 
 [z1,z2] = InputSegment("z1", "z2");
 axis([x1,x2,y1,y2,z1,z2]);
 replot;
endif;

endfunction;


function TitlePlot(n)

header = input("Title = ", "s");
title(header);

xl     = input("Label for x axis = ", "s");
xlabel(xl);

yl     = input("Label for y axis = ", "s");
ylabel(yl);

if (n >= 3)
 zl    = input("Label for z axis = ", "s");
 zlabel(zl);
endif

replot;

endfunction


function LoadDefinition()
input("Type file name (without the extension .m).");
endfunction 


disp(" INTRODUCTION"),
disp("To have something to plot, write a file with your Octave definitions of ")
disp("functions (extension .m). Look through Octave definitions of functions ")
disp("in this file to learn how to write down your file.")
disp("  *A few preloaded examples available right now. Check in items 1-5.*  ")
disp(" ")

% "n" marks the "New Plot" phase.
m = "n";

while (m == "n")

disp(" MAIN MENU"),
disp("0.Load a definition file."), 
disp("1.Draw a line    y       = Y(x)."),
disp("2.Draw a line    (x,y)   = F(t)"),
disp("3.Draw a line    (x,y,z) = F(t)."),
disp("4.Draw a surface z       = Z(x,y)."),
disp("5.Draw a surface (x,y,z) = F(u,v)."),
disp("6.Quit.")
disp(" ")
n = input("Select item number n = ");

% We set axis to autoscaling
axis();

if (n == 0) 
 LoadDefinition(); 		
endif;
if (n == 6) 
 m = "q";			
endif;

if ((1 <= n) & (n <= 5))
if (n == 1)
 disp ("Plot of y = Y(x) over [x1,x2], using n subdivision."),  
 Y = input("Preloaded example Y=Cube. Insert the name of Y= ","s");
endif;
if (n == 2)
 disp ("Plot over [t1,t2], using n subdivisions, of (x,y) = F(t)."),  
 disp("Preloaded example F = Circle.")    
 F = input("Insert the name of F = ", "s");
endif;
if (n == 3)
 disp ("Plot over [t1,t2], with n subdivisions, of (x,y,z) = F(t)."),  
 disp("Preloaded example F=Elicoid.")
 F = input("Insert the name of F= ", "s");
endif;
if (n == 4)
 disp ("Plot over [x1,x2]X[y1,y2], with nxn subdivisions, of z = Z(u,v)."),  
 Z = input("Preloaded example Z=Wave1. Insert the name of Z= ", "s");
endif;
if (n == 5)
 disp ("Plot over [u1,u2]X[v1,v2], with nxm subdivisions, of (x,y,z)=F(u,v)."),  
 disp("Preloaded example F=Sphere.")
 F = input("Insert the name of F = ", "s");
endif;

% "i" marks the "Independent variable setting" phase.
m = "i";

while (m == "i")

if (n == 1) LinePlot(Y);		endif;
if (n == 2) LineParPlot(F);		endif;
if (n == 3) Line3DPlot(F);		endif;
if (n == 4) SurfacePlot(Z);		endif;
if (n == 5) SurfaceParPlot(F);		endif;
if (n == 6) m = "q";			endif;

 % "d", "t", "s" mark a short phase of fine tuning of the same plot. 
 m = "d";
 while ((m == "d") | (m == "t") | (m == "s"))
  disp(" PLOT MENU")
  disp("D (or I) = change dependent (or independent) variable range.")
  disp("T= add/change title and labels. S=Save, N=new plot, Q=quit.")
  m = input(" Choice =  ","s");
  m = tolower(m);
  if (m == "d") ChangePlot(n); 	endif;  
  if (m == "t") TitlePlot(n);	endif;
  if (m == "s") SavePlot(n);	endif;
 endwhile;
endwhile;
endif;
endwhile