clc;clear;close all;
% Sphere Project
 
% Z axis
dz=10; nz=500; oz=10;
 
nx=256;
ny=256;
a1max=65;
a2max=35;
LL=linspace(-a1max,a1max,nx)/180*pi;
PP=linspace(-a2max,a2max,ny)/180*pi;
ZZ=[oz:dz:(nz-1)*dz+oz];
dx=oz*(LL(2)-LL(1));
dy=oz*(PP(2)-PP(1));
ox=oz*LL(1);
oy=oz*PP(1);
 
% . .  Set up coords
phi0=acos(2/pi);
 
for jj=1:ny
    disp(['Working on ',num2str(jj),' of ',num2str(ny)])
    for ii=1:nx
        AA=acos( cos(PP(jj))*cos(LL(ii)/2));
        if (sin(AA) == 0)
            WW=0;
        else
            WW = 1/sin(AA);
        end 
        tmp1 = (2*LL(ii)/pi + 2*WW*AA*cos(PP(jj))*sin(LL(ii)/2))/2;
        tmp2 = ( PP(jj)+2*WW*AA*sin(PP(jj)) )/2;
        for kk=1:nz
            xx(kk,ii,jj)=ZZ(kk)*tmp1;
            yy(kk,ii,jj)=ZZ(kk)*tmp2;
            if (1-tmp1^2-tmp2^2 > 0) 
                zz(kk,ii,jj)=ZZ(kk)*real(sqrt(1-tmp1^2-tmp2^2));
            else
                zz(kk,ii,jj)=-ZZ(kk)*real(sqrt(tmp1^2+tmp2^2-1));  
            end
        end
    end
end
 
 disp('PRINTING FILE')
ttest=1;
if (ttest==1) 
    figure(3)
    surf(squeeze(yy(20,:,:)),squeeze(xx(20,:,:)),-squeeze(zz(20,:,:)));
    alpha(0.75);
    axis equal;
    xlabel('Inline Direction (m)');
    ylabel('Crossline Direction (m)');
    zlabel('Depth (m)');
    view(45,20)
   % print -depsc2 Starting.ps
end 
 
rays = cat(4,cat(4,xx,yy),zz);

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % . . Coordinate system
  eval(['D1.history_file=''WTRays.HH'' '])
  D1.n=[nz nx ny 3];   D1.o=[oz ox oy 0];
  D1.d=[dz dx dy 1];   D1.esize=4;
  eval(['D1.in=''/net/koko/jeff/WTRays@'' ']);
  D1.data_format='xdr_float';
  D1.data=rays;
  eval(['sepwrite(D1,''WTRays.HH'') '])
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

