# Sample parametric plot

# You can run this by typing at the prompt: octave -q pplot.m

outfile = "pplot.dat";
global z = 1.3;
a_vec = linspace(-pi/2, pi/2, 100);

function r = g(x)
  global z;
  r = 2*sin(x)./sqrt((z + sin(x)));
endfunction

function r = G(a)
  r = quad("g", a, pi/2);
endfunction

function r = h(x)
  global z;
  r = 2*sin(x).*sqrt((z + sin(x)));
endfunction

function r = H(a)
  r = quad("h", a, pi/2);
endfunction

points = [];
for a = a_vec
  points = [points G(a)/sqrt(z + sin(a)) H(a)/(z + sin(a))**(3/2)];
endfor

fd = fopen(outfile, "wt");
fprintf (fd, "%14.6f %14.6f\n", points);
fclose(fd);
disp("Done!");