function meshGWPulsar(n,lambda,a,b,rMin,rMax,tOverT)
% usage: (n,lambda,rMin,rMax,tOverT)
% Generates surf plot of GWs from rotating elliptically shaped gravity-well (i.e., a GW pulsar).
% n = number of grid points for n X n grid
% lambda = wavelength of the waves
% tOverT = time for which to make the mesh plot divided by period of GWs
% a,b = constants in equation for rotating ellipse, xrot^2/a^2 + yrot^2/b^2 = 1
% rMin = minimum radius to show with waves
% rMax = maximum radius to show with waves

if (a > rMin || b > rMin)
    disp('a and b must be <= rMin')
    return;
end

c = 1.0; % Think of time in years, distance in light-years
k = 2*pi/lambda;
w = k*c;
T = 2*pi/w;
t = tOverT*T;
Omegat = (w/2.0)*t; % angular velocity of rotating ellipse at half GW freq, times t.
a2 = a*a;
b2 = b*b;
delta = rMax/n;
xVec = -1.0*rMax + (0:2*n)*delta;
yVec = -1.0*rMax + (0:2*n)*delta;
for i = 1:2*n
    for j = 1:2*n
        x = xVec(j); % x value for column j
	y = yVec(i); % y value for row i
	r = sqrt(x*x + y*y);
	phi = atan2(y,x);
	if r >= rMin
	  z(i,j) = (1.0/r)*cos(w*t - k*r - 2*phi);
	else
	  %rotate x and y by Omegat:
	  xrot = x*cos(Omegat) + y*sin(Omegat);
	  yrot = -1.0*x*sin(Omegat) + y*cos(Omegat);
	  d2 = (xrot*xrot/a2) + (yrot*yrot/b2);
          if (d2) <= 1
 	    z(i,j) = (d2 - 1.0)/rMin; % inside the rotated ellipse draw the gravity well of the star
	  else
 	    z(i,j) = 0.0;
	  end 
	end
    end
end
surf(z)
axis off
view([0 0.5*rMax 2*rMax]);
shading interp
%hidden off
