% contour plot of V(x,y)=1/2 y^2 - kappa*cos(x)
%
% jbr, 9/3/2011
%
n = 101;
x = linspace(-5, 5, n);
y = linspace(-5, 5, n);
kappa = 2;
xx = kron(x, ones(n,1));
yy = kron(y', ones(1,n));
yy2 = yy.*yy;
z = 0.5*(yy2) - kappa*cos(xx);
val = [-1, 0, 1, 2, 3, 5, 8, 12];
[c, lev] = contour(x,y,z,val);
axis([-5, 5, -5, 5])
xlabel ("$q$");
% avoid rotating y-axis label
text (-5.75, 0, "$p$")
% write contour data for gnuplot
col = 1;
myfile = fopen('heteroclinic.dat', 'w');
while col <= size(c, 2)
coln = col + c(2,col);
data = c(:,col+1:coln);
% print this contour data and move to next contour line
fprintf(myfile, '%e %e \n', data)
fprintf(myfile, '\n\n')
col = coln + 1;
end
fclose(myfile);