function lmcosip=plm2pot(lmcosi,r,GM,a,wat)
% lmcosip=PLM2POT(lmcosi,r,GM,a,wat)
%
% Converts the real spherical harmonic coefficients expressing a gravity
% field to actual values of the gravitational potential. If the (0,0)
% term is missing, can add in the reference constant for that as well, and
% also supplies the zeros for (1,0) and (1,1).
%
% INPUT:
%
% lmcosi        [l m Ccos Csin] degrees, order, coefficients
% r             Requested radius, in m
% GM            Gravitational coefficient times mass reference constant
% a             Reference radius
% wat           1 gravitational potential
%               2 free-air gravity anomaly
%               3 geoid anomaly
%
% OUPUT:
%
% lmcosip       [l m Ccos Csin] degrees, order, coefficients of output
%
% Last modified by fjsimons-at-alum.mit.edu, Feb 13th 2004

% See Research Book V page 84, Essentials 2001 and notes with egm96_?

defval('GM',fralmanac('GM_EGM96','Earth'))
defval('a',fralmanac('a_EGM96','Earth'))
defval('r',a)
defval('wat',1)

GMr=GM/r;

% Blakely Eq (7.2)
el=lmcosi(:,1);
arl=(a/r).^el;

lmcosip=lmcosi;

switch wat
 case 1 % Potential - Blakely (7.2)
  fact=GMr*[arl arl];
  disp('Calculating gravitational potential')
 case 2 % Free air gravity anomaly, NOT the gravity disturbance!!!
  fact=GMr*[arl arl].*[el-1 el-1]/r;
  disp('Calculating free-air gravity anomalies')
 case 3 % Geoid anomaly (undulation)
  fact=a;
end

lmcosip(:,3:4)=fact.*lmcosi(:,3:4);
%lmcosip=taper(lmcosip,50,60,'cosine');
%disp('Do not forget this is tapered')
%load taperprop
%lmcosip=taper(lmcosip,tpr(n,1),tpr(n,2),'cosine');

% Mean and center of mass were zero
if lmcosi(1)==2
  lmcosip=[0 0 GMr 0 repmat(0,1,size(lmcosi,2)-4);
	   1 0 0   0 repmat(0,1,size(lmcosi,2)-4);
	   1 1 0   0 repmat(0,1,size(lmcosi,2)-4);
	   lmcosip];
elseif lmcosi(1)==1 & lmcosi(1,2)==0
  lmcosip=[0 0 GMr 0 repmat(0,1,size(lmcosi,2)-4);
	   lmcosip];
end

% Now this is complete. Try to take out the C(2,0) terms
% Take out the C(0,0) term again too, to demean
lmcosip(1,3)=0;
lmcosip(4,3)=0;