%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [rlat,rlon,rlonc] = radii(lat)

R = 6378137;
r = 6356752.3142;
e = 0.0818191908426;
f = 1/298.257223563;

esq = e*e;
Rsq = R*R;
rsq = r*r;
latr = deg2rad(lat);
slat = sin(latr);
clat = cos(latr);
tlat = tan(latr);
sinsqlat = slat.*slat;
cossqlat = clat.*clat;
tansqlat = tlat.*tlat;
tmp1 = sqrt(1-esq*sinsqlat);
rlat = R*(1-esq)./tmp1.^3;
rlon = R./tmp1;
rlonc = rlon.*clat;