function [xx,yy]=projocta(lono,lato,R,nsect);
% function [xx,yy]=projocta(lono,lato,R,nsect);
% b klinger, 7/15/12
% Project (longitude,latitude)(lono,lato) on to (x,y) 
% 2D arrays (xx,yy) in order to create "OctaGlobe" polyhedron with 
% longitude and latitude profiles that are octagons and
% vertices on a sphere of a given radius R (defaults for no R or R==0).
% nsect=[N S] sectors to place polar octagons 
% (range 1 to 8, default=[8 4])

%          MODIFICATIONS
% 8/21/12 nsect option
% 5/1/16 create default value of R
% 5/3/16 can take either 1D or 2D (lono,lato)
% 5/7/16 nsect=[0 0] allows nsect to take default values

% Parameters and arrays based on (lono,lato).

if min(size(lono))==1 & min(size(lato))==1
   [lon,lat]=meshgrid(lono,lato);
 else
   lon=lono;
   lat=lato;
end
nnx=size(lon);
nny=size(lat);
if nnx~=nny
   disp('arrays have inconsistent sizes')
   return
end
nx=nnx(2);
ny=nnx(1);
xx=zeros(ny,nx);
yy=zeros(ny,nx);
lon0=22.5+45*floor(lon/45);
dlon=mod(lon,45)-22.5;
costh=cos(abs(dlon)*pi/180);

% Choose orientation and location of polar octagons
% (There are 8 tropical/midlatitude sectors, starting at
% longitude 0-45E; nsect tells which of these to attach
% Arctic and Antarctic octagons.)

%% angN=3.5*45; xsectN=1.5;   N sample 1
%% angS= .5*45; xsectS=5.5;   S sample 1
%% angN=1.5*45; xsectN=7.5;   N sample 2
%% angS=2.5*45; xsectS=3.5;   S sample 2

if ~exist('R')
   R=1 ./sqrt(sqrt(2)*(sqrt(2)-1));   % choose R to make A=1
end
if R==0
   R=1 ./sqrt(sqrt(2)*(sqrt(2)-1));   % choose R to make A=1
end
if ~exist('nsect')
   nsect=[8 4];
end
if nsect==[0 0]
   nsect=[8 4];
end
angN=45*(nsect(1)+1.5);
xsectN=(nsect(1)-.5);
angS=45*(6.5-nsect(2));
xsectS=(nsect(2)-.5);

% cos/sin arrays in order to 
% attach Arctic face to NH midlat face 2 (45-90 deg) and
% Antarctic face to SH midlat face 6 (225-270 deg)

coslonN=cos(pi*(lon-angN)/180);
sinlonN=sin(pi*(lon-angN)/180);
coslonS=cos(pi*(-lon-angS)/180);
sinlonS=sin(pi*(-lon-angS)/180);

% Size parameters for tropical faces (A,B), mid-latitude
% faces (C,D), and polar faces (E).

A=sqrt(sqrt(2)*(sqrt(2)-1))*R;
B=R/sqrt(2);
C=sqrt(5/4-1/sqrt(2))*R;
D=sqrt(3/2-sqrt(2))*R;
E=R/sqrt(2);

% Loop through all latitudes and construct (x,y) arrays

ijT=find(abs(lat)<=22.5);
xx(ijT)=B*lon(ijT)/45;
yy(ijT)=A*(lat(ijT)/45);

ijM=find(abs(lat)>22.5 & abs(lat)<=67.5);
dlatr=(abs(lat)-22.5)/45;
yy(ijM)=sign(lat(ijM)).*((A/2)+C*dlatr(ijM));
xx(ijM)=B*lon0(ijM)/45+(dlon(ijM)/45).*(B-(B-D)*dlatr(ijM));

ijA=find(lat>67.5);
r=E*((90-abs(lat))/45)./costh;
xx(ijA)=xsectN*B+r(ijA).*coslonN(ijA);
yy(ijA)=(A/2+C+E/2)+r(ijA).*sinlonN(ijA);

ija=find(lat<-67.5);
r=E*((90-abs(lat))/45)./costh;
xx(ija)=xsectS*B+r(ija).*coslonS(ija);
yy(ija)=-(A/2+C+E/2)+r(ija).*sinlonS(ija);

