function [w, g] = rotpro(nl,nb) 
% Phase-dependent rotation profile.
% Needs lots of resolution in longitude (L) direction.
%
% D. Holmgren holmgren@phobos.astro.uwo.ca

incl = input('Enter i: ');
phase = input('Enter phase: ');
nw = input('Number of lambda bins: ');
x = ones(nl*nb,1);
phase = phase * 2*pi;
ci = cos(incl*pi/180);
si = sin(incl*pi/180);
L = linspace(0,2*pi,nl)'; dl = 2*pi/nl;
B = linspace(-pi/2,pi/2,nb)'; db = pi/nb;
g = zeros(nw,1); w = linspace(-1,1,nw)';
u = 0.6;
for m = 1:nb
sb = sin(B(m));
cb = cos(B(m));
ds = db*dl*cb;
 for l = 1:nl
  j = m + nb*(l-1);
  cth = sb*ci + cb*si*cos(L(l) + phase);
  brk = 1;
  if cth <= 0, brk = 0;, end
  y = cb * sin(L(l) + phase);
  A = (1. - u + u*cth)*cth;
  xa = x(j)*A*brk*ds;
%  khi = find(y >= w); kw = max(khi);
  for k = 1:nw-1
  if y < w(k+1) & y >= w(k),
   kw = k;
  end
  end
  g(kw) = g(kw) + xa;
 end
end
g = g ./ sum(g);