% Plot contours in SPEM orthogonal curvilinear grid coordinates
% Test m-file
% John Wilkin 12 Apr 95

[cs,ha]=contour(h,30);

for i=1:length(ha),
  x=get(ha(i),'xdata'); y=get(ha(i),'ydata'); 

  % x is i index of spem coords
  % y is j index of spem coords
  
  [m n] = size(rholon);
  X = interp2(1:n,1:m,rholon,x,y);
  Y = interp2(1:n,1:m,rholat,x,y);
  
  % X is longitude at (i,j)
  % Y is latitude at (i,j)

  set(ha(i),'xdata',X); set(ha(i),'ydata',Y);
end
axis auto

