%  this MATLAB file makes contour and 3-d plots

%  of a user-specified function f(x,y)

% set and b
a=2;b=2;

% set Vo
Vo=1;

%set the number of grid points
npts=50;

%set minimum x
xmin=-b;

%set maximum x
xmax = b;

%set minimum y
ymin=0;

%set maximum y
ymax = a;

%set delta x
dx = (xmax-xmin)/npts;

%set delta y
dy = (ymax-ymin)/npts;

%x-range of plotting points
x = xmin:dx:xmax;

%y-range of plotting points
y = ymin:dy:ymax;

%Enter the maximum n to which to take the sum
nmax=input(' n max - ')

%loop over the grid and evaluate f(x,y)

nn= length(x) ;
mm= length(y) ;
clear V ;

% store the functions that will be used over and over
% so that the code will run fast

% fy is a matrix containing sin(ky) for all of
% the k-values and for all of the y's on the grid

% fx is a matrix containing cosh(kx)/cosh(ka) for all of
% the k-values and for all of the x's on the grid

fx=zeros(nmax,nn);
fy=zeros(nmax,mm);

for n=(1:2:nmax)
   k=n*pi/a;
   fx(n,:) = cosh(k*x)/cosh(k*b);
   fy(n,:) = sin(k*y);
end

% build the potential V(x,y) at the grid points
% i moves across the grid in x

for i=1:nn
   % j moves across the grid in y
   for j=1:mm
      V(i,j)=0;

      % do the terms in the sum
      for n=1:2:nmax
         V(i,j) = V(i,j) + fx(n,i)*fy(n,j)/n;
      end
   end
end

%  multiply by the physical constant (top of page 134)
V = 4/pi*Vo*V;

% look at V in the middle
V((nn+1)/2,(mm+1)/2)

%make a surface plot of the potential
 mesh(x,y,V')
 xlabel('x');ylabel('y');title('V(x,y)')

% viewing altitude and azimuth in degrees
alt=50;azi=-30;
view(azi,alt)

pause;

%make a 20-contour plot
contour(x,y,V',40)
 xlabel('x');ylabel('y');title('V(x,y)')