**************************************************** ;

goptions reset=all vsize=7.1 in  ;
* filename grafout 'c:grafout\globe5.ps';
 ** goptions reset to file just before proc greplay ***;

goptions targetdevice=pslepsf;
proc greplay nofs igout=work.gseg; delete _all_;

options nomacrogen ;  goptions lfactor=2 aspect=0;
data twoax; input x y; cards;
 -.1  -.2          These will delimit the 2-D plots
  1.1   1.2
 ;
data globe;
 input id $ x y z; nnn+1;
  if y=99 then y=sqrt(1-x*x);
  if x=99 then x=sqrt(1-y*y);
  put 'X=' x ' Y=' y;
cards;
Origin  0   0   0
 X1     1   0   0      These points will provide the ends
 X2     0   1   0       of the axes in our plot.
XX1     .85  99 0
XX2      99 .85 0
top     0   0   .2
bot .00001  0  -.2     <----- Jitter
 ;
data ball;   length color $8. function $5. id $10. style $12.;
 style = 'duplex ';
 pi = 4*atan(1); put 'Using pi = ' pi;
 a=.85; b=.53; rmax=.12;  drop a b;
 xsys='2'; ysys='2'; zsys='2';
   when='a';
                    /*---------------------------------------------
                    |  Make the ball and its equatorial circle     |
                     ---------------------------------------------*/

  ****** Slices on X are circles *******;
               id = 'ball   ';                     color = 'cyan';

  do x = -.65,-.3,0,.3, .65;  r=sqrt(1-x*x);
  function = 'move';
  do arg=0 to 2 by .1;
  y = r*cos(arg*pi); z=r*sin(arg*pi);
  output; function = 'draw'; end; end;

  ***** Slices on Y are circles ******;
  do y = -.65,-.3,0,.3 ,.65;  r=sqrt(1-y*y);
  function = 'move';
  do arg=0 to 2 by .1;
  x = r*cos(arg*pi); z=r*sin(arg*pi);
  output; function = 'draw'; end; end;

  ***** Slices on Z are circles too ****;
  do z = -.9 to .9 by .1;
  r=sqrt(1-z*z);
  function = 'move';
  do arg=0 to 2 by .1;
  y = r*cos(arg*pi); x=r*sin(arg*pi);
  output; function = 'draw'; end; end;

  **** Projected onto floor (Z=0) ***;
          id='circle';
  z=0; function = 'move';
  do arg= 0 to 2 by .01;                     color = 'blue';
  y =  cos(arg*pi); x= sin(arg*pi);
  output; function = 'draw'; end;





                   /*-------------------------------------------
                   |  Make the X and W axes and their labels    |
                    -------------------------------------------*/

  **** X-space axes orthogonal *****;
          id = 'X_Axes';
  z=0; x=0; y=0; function='move';            color='green';
  output;
  x=1; y=0; function = 'draw'; output;
  x=0; y=0; function = 'move'; output;
  x=0; y=1; function = 'draw'; output;
  x=0; y=0; function = 'move'; output;
  z=-.2; function = 'draw'; output; z=.2; output;  z=0; output;
  ** Axis labels **;
  x=1.04; y=0;  function='label'; position='5'; text = 'X ';
       size=2.5; output;
       x=x+.02; text='1'; position='9'; size=1.5; output;
  x=-.02; y=.95 ; function='label'; position='5'; text = 'X ';
       size=2.5; z=-.05; output;
       x=x+.02; text='2'; position='9'; size=1.5 ;  output;
  function='draw'; position = '5'; z=0;

  **** X-space axes not orthogonal *****;
                  id='W_Axes';
  line=1;                                   color='gray50';
  x=0; y=0; output; x= a; y= b  ; output;
  x=0; y=0; output; x= b; y= a  ; output;
  ** Axis labels **;

  x=a*1.1; y=b*1.1; function='label'; text = 'W'; size=2.5; output;
     x=x+.02;
     position='9'; text='1'; size=1.5 ; output; position = '5';
  x=b*1.1; y=a*1.1; z=z+.04; function='label'; text = 'W'; size=2.5; output;
     x=x+.02;
     position='9'; text='2'; size=1.5 ; output; position = '5';
  z=0;
  line=1;

                    /*------------------
                    |  Make the fences  |
                     ------------------*/


  **** Project some circle points onto these axes ***;
             id = 'Fence';
  function = 'move'; z=0;

  %macro fence(a,b,c1,ht,line ) ;

     color= " &c1 "; line=&line; function = "move";
     *** The point (x0,y0) is on the circle ****;
          z=0; x=x0; y=y0; output;

     *** Draw to axis whose slope is a/b  *******;
  x =  &a*(&b*y0 - &a*x0)/(&b*&b-&a*&a);
  if x>0 then y =  &b*x/&a; function ="draw";  output;
  x0a = x; y0a = y; **<----- axis projection = end of fence ****;


     *** The point (x1,y1) is on the circle ****;
           x=x1; y=y1; function="move";  output;

     *** Draw to the  axis *******;
  x =  &a*(&b*y1 - &a*x1)/(&b*&b-&a*&a);
  if x>0 then y =  &b*x/&a;  function ="draw";  output;
     *** Make the fence ****;
  z = &ht; idg=1;  output;
    x=x0a; y=y0a; idg=2; output;
  z=0; idg=3;  output; idg=.;

  %mend fence;


   x0=.5+rmax*cos(3.6*pi/6); y0=.5+rmax*sin(3.6*pi/6);
   x1=.5+rmax*cos(-.5*pi/6); y1=.5+rmax*sin(-.5*pi/6);

 id='Wfence1';  %fence( a,b, pink, -.115,8)
 id='Xfence1';  %fence( 1,0, brown,-.08,1)
 id='Wfence2';  %fence( b,a, pink,  .115,8)
 id='Xfence2';  %fence( 0,1, brown, .08,1)


** Macro for letter with subscript and hat **;
  %macro greek(x,nn);
  if z>0 then fact= 1.85; else fact=1.45;
  z = z*fact; function='label';
  position = '5'; style = 'cgreek'; text = &x ; size=2.5; output;
  position = 'B'; style='special'; text = '1';           output;
  position = '9'; style = 'cgreek'; text = &nn; size=1.5  ; output;
  z=z/fact; position = '5';
  function = 'draw';
  %mend greek;

**** Center the ball at right point ******;
**** Add the Greek letters with subs *****;
data anno; length style $ 8.;  set ball;
  if id='ball' or id = 'circle' then do;
      x=rmax*x+.5;
      y=rmax*y+.5;
     if id = 'ball' then z = rmax*z  /*+.5*/ ;
  end; output;
  if id='Xfence1' and idg=1 then do; %greek('b','11'); end;
  if id='Xfence1' and idg=2 then do; %greek('b','12'); end;
  if id='Xfence2' and idg=1 then do; y=y-.03; %greek('b','21');y=y+.03; end;
  if id='Xfence2' and idg=2 then do; y=y+.03; %greek('b','22');y=y-.03; end;
  if id='Wfence1' and idg=1 then do; %greek('g','11'); end;
  if id='Wfence1' and idg=2 then do; %greek('g','12'); end;
  if id='Wfence2' and idg=1 then do; z=z+.03; x=x-.01;
                                     %greek('g','21');z=z-.03; x=x+.01; end;
  if id='Wfence2' and idg=2 then do; %greek('g','22'); end;

proc freq data=anno; tables id*color;


data anno1; set anno;  drop line;
if color = 'cyan' then color = 'gray50';
else  color = 'black';
  if id='Xfence2'  and function = 'label' then z = z+.02;
  if id='Wfence1'  and function = 'label' and idg= 1 then do;
     x = x-.025 ; z=z-.01; end;
  if id='Wfence2'  and function = 'label' and idg= 2 then do;
     x = x-.04 ;   end;
  if id='Wfence2'  and function = 'label' and idg= 1 then do;
     z = z-.04 ;   end;


proc g3d data=globe annotate=anno1;
  scatter x*y=z/shape = 'balloon' size=.2 color='black'
  rotate = 20 tilt = 50 noneedle noaxes nolabel;
  title  h=3 j=center f=duplex '(a)';
proc g3d data=globe annotate=anno1;
  scatter x*y=z/shape = 'balloon' size=.2  color='black'
  rotate = 5 tilt =  5 noneedle noaxes nolabel;
proc g3d data=globe annotate=anno1;
  scatter x*y=z/shape = 'balloon' size=.2  color='black'
  rotate = 75 tilt =  90 noneedle noaxes nolabel;


data anno2; set anno;
  if id ne 'W_Axes' and id ne 'ball  ' and substr(id,1,4) ne 'Wfen';
  if id = 'X_Axes' and function= 'label' then y=y*1.3;
  if id='Xfence1'  and function = 'label' then do;
     y = y-.07; end;
  if id='Xfence2'  and function = 'label' then do;
     x = x-.07 ; y=y+.05; end;
proc freq; tables color*id;  title 'Anno2 data';
proc print data=anno2;
  var function id x y z color;
  where id ne 'ball';

data anno3; set anno;
  if color ne 'cyan' and color ne ' brown ' and color ne 'green';
  if x=0 and y = 0 then function = 'move';
  if id='Wfence1' and function = 'label' then y = y - .09 ;
  if id='Wfence2' and function = 'label' then do;
            y = y + .1; x=x-.07; end;

  proc freq; tables color*id; title 'Anno3 data';
data anno2; set anno2; color = 'black'; drop line;
data anno3; set anno3; color = 'black'; drop line;


                 goptions aspect=1.2;

data anno2; set anno2;
 if y>1.234 then do; put @20 x y;  y = 1; x=x+.05; output; end;
 else output;

proc gplot data=twoax annotate=anno2;
   plot y*x/ noaxes nolegend ;
 symbol1 v=none i=none c=purple;
 title  j=center h=2.6 f=duplex '(b)';
proc gplot data=twoax annotate=anno3;
   plot y*x/ noaxes nolegend  ;
 symbol1 v=none i=none c=purple;
 title  j=center h=2.6 f=duplex '(c)';

  /*
goptions
 dev=pslepsf gsfmode=replace gsfname=grafout targetdevice=pslepsf;
           */
 footnote ' ';

  proc greplay tc=tempcat nofs igout=work.gseg;
  tdef globe des='Globe'

  1/ llx = 15   lly = 55
     ulx = 15   uly = 96
     urx = 85   ury = 96
     lrx = 85   lry = 55


   2/llx =  5   lly =  1
     ulx =  5   uly = 60
     urx = 50   ury = 60
     lrx = 50   lry =  1

   3/llx = 50   lly =  1
     ulx = 50   uly = 60
     urx = 95   ury = 60
     lrx = 95   lry =  1


;

   template = globe;

  treplay 1:g3d  2:gplot 3:gplot1;
  replay template;

run ;