%!
%----------------------------
%-----define constants-------
%----------------------------
/width              300   def  %screen
/height             300   def  %

/sphere_x             0.0 def  %sphere object
/sphere_y             0.0 def  %
/sphere_z           200.0 def  %
/sphere_radius      150.0 def  %
/sphere_ambient       0.2 def  %
/sphere_specular      0.4 def  %
/sphere_specular_n   10.0 def  %
/sphere_diffuse       0.4 def  %

/light_vx -1.0 1.732  div def  %light vector 
/light_vy  1.0 1.732  div def  %
/light_vz  1.0 1.732  div def  %

/eye_x                0.0 def  %place of eye
/eye_y                0.0 def  %
/eye_z          200.0 neg def  %

%----------------------------
%-----define procedure-------
%----------------------------
/Abs
  {
    dup
    0 ge
    {
    }
    {
      neg    
    } ifelse
  } def

/Root 
% stack top become from a^2 to a  
  {
    dup
    20                         % 20 times loop
    {
      /Root_backup_x exch def  % Root_backup_x = x
      dup                      % a a
      Root_backup_x            % a a x
      2.0 div                  % a a (x/2.0)
      exch                     % a (x/2.0) a 
      Root_backup_x 2.0 mul    % a (x/2.0) a (x*2.0)
      div                      % a (x/2.0) a/(x*2.0)
      add                      % a (x/2.0)+a/(x*2.0)
    } repeat
    exch
    pop
  } def

/Raytrace
  {   
    /i height 2.0 div neg def    
    height                     % height times loop
    {
      /j width 2.0 div neg def   
      width                    %width times loop
      {
        /E_x   j eye_x sub def
        /E_y   i eye_y sub def
        /E_z 0.0 eye_z sub def
        /E_len E_x dup mul     
               E_y dup mul     
               add             
               E_z dup mul     
               add             
               Root           
               def             
        /E_x E_x E_len div def   
        /E_y E_y E_len div def
        /E_z E_z E_len div def
        /A 1 def    
        /B E_x eye_x sphere_x sub mul 
           E_y eye_y sphere_y sub mul
           add
           E_z eye_z sphere_z sub mul 
           add
           def
        /C eye_x sphere_x sub dup mul
           eye_y sphere_y sub dup mul
           add
           eye_z sphere_z sub dup mul
           add
           sphere_radius dup mul 
           sub
           def
        /D B dup mul 
           A C mul
           sub
           def
        gsave
        D 0.0 ge        % if (D >= 0.0) 
        {               % then       
           /t B neg D Root sub def
           t 0.0 ge
           {
              /now_x eye_x E_x t mul add def
              /now_y eye_y E_y t mul add def
              /now_z eye_z E_z t mul add def
              /n_vx now_x sphere_x sub def
              /n_vy now_y sphere_y sub def
              /n_vz now_z sphere_z sub def
              /n_len n_vx dup mul
                     n_vy dup mul
                     add
                     n_vz dup mul
                     add
                     Root
                     def
              /n_vx n_vx n_len div def
              /n_vy n_vy n_len div def
              /n_vz n_vz n_len div def
              /nl   n_vx light_vx neg mul
                    n_vy light_vy neg mul
                    add
                    n_vz light_vz neg mul
                    add
                    def
              /r_vx light_vx nl Abs div n_vx 2.0 mul add def
              /r_vy light_vy nl Abs div n_vy 2.0 mul add def
              /r_vz light_vz nl Abs div n_vz 2.0 mul add def
              /r_len r_vx dup mul 
                     r_vy dup mul
                     add
                     r_vz dup mul
                     add
                     Root
                     def
              /r_vx r_vx r_len div def
              /r_vy r_vy r_len div def
              /r_vz r_vz r_len div def
              /cos_gamma E_x neg r_vx mul
                         E_y neg r_vy mul
                         add
                         E_z neg r_vz mul
                         add
                         def
 	      /result sphere_ambient 1.0 mul
                      sphere_diffuse 1.0 mul nl mul
                      add
                      def
              cos_gamma 0 gt
              {
                 /result result
                         1.0 sphere_specular mul
                         10
                         {
                            cos_gamma mul
                         } repeat
                         add
                         def
              }
              {
              } ifelse
              result setgray
           }
           {
              0.0 setgray
           } ifelse
        }
        {               % else
           0.0 setgray
        } ifelse
        1 0 rlineto
        0 1 rlineto
        1 neg 0 rlineto
        closepath
        fill
        grestore        
        1 0 rmoveto
        /j j 1 add def
      } repeat
      width neg 1 neg rmoveto
      /i i 1 add def           
    } repeat
    stroke
  } def

%----------------------------
%-----main program-----------
%----------------------------
newpath
150 600 moveto
Raytrace
showpage
quit