#!/usr/bin/env python

from math import *
M_PI = pi
true=1
false=0

from epix import *

import sys

####################################################################
# Converted from clipping.xp
####################################################################

# Style parameters
# camera location (in spherical coordinates); must be in first orthant
VIEWPT=sph(4, M_PI/6, M_PI/6)

LIGHT=-E_1-E_2-2*E_3 # location of light, for shading

def color_coord():
  rgb(0.7, 0.9, 0.5)
def color_axis():
  rgb(0.8, 0.2, 0.9)
def color_xslice():
  red()
def color_yslice():
  blue()
def ps_axis_color():
  print "\n\\newrgbcolor{axis_color}{0.8 0.2 0.9}%"

MESH=8 # number of coordinate grid squares

# location of tangency point
x_0=0.625
y_0=0.5

z_0=0.25 # height of top of slicing planes

# function to be graphed
def f(x,y):
  return 0.75*y*(y*y-3*x*x)

# Auxiliary internal parameters
N=2*MESH # number of surface mesh subdivisions
MAX=1 # maximum coordinate

# constant y slice of the graph of f
def Df_x(t):
  return P(t+x_0, y_0, f(t+x_0, y_0))

# constant x slice of the graph of f
def Df_y(t):
  return P(x_0, t+y_0, f(x_0, t+y_0))

# shading density function and helpers
EPS=0.001 # for approximating derivatives
def g(i):
  return MAX*i*1.0/N # int -> grid coordinate

# 0.5*(1-cos^2(t)), t=angle from camera to surface normal
def density(i, j):
  # compute unit surface normal
  f_x = 0.5*(f(g(i)+EPS, g(j)) - f(g(i)-EPS, g(j)))/EPS
  f_y = 0.5*(f(g(i), g(j)+EPS) - f(g(i), g(j)-EPS))/EPS

  normal = P(-f_x, -f_y, 1)
  normal *= 1/sqrt(normal|normal)

  # convert cos^2 to a gray density
  #  return  0.5*(1-pow(normal|VIEWPT, 2)/(VIEWPT|VIEWPT))
  return  0.625*(1-pow(normal|LIGHT, 2)/(LIGHT|LIGHT))

# mesh location -> draw surface element
def draw_mesh_quad(i, j):
  gray(density(i,j))
  quad(P(g(i),   g(j),   f(g(i),   g(j))), 
       P(g(i+1), g(j),   f(g(i+1), g(j))),
       P(g(i+1), g(j+1), f(g(i+1), g(j+1))), 
       P(g(i),   g(j+1), f(g(i),   g(j+1))))

bounding_box(P(-2,-2), P(2,1))
picture(2,1.5)
unitlength("1.25in")

use_pstricks()
begin()
viewpoint(VIEWPT)
use_pstricks(false)

clip_box(P(-2, -2, -1), P(2, 2, 0.75))

# coordinate grids
color_coord()
grid(P(-MAX,-MAX,-MAX), P(-MAX, MAX, MAX), MESH, MESH)
grid(P(-MAX,-MAX,-MAX), P( MAX,-MAX, MAX), MESH, MESH)
grid(P(-MAX,-MAX,-MAX), P( MAX, MAX,-MAX), MESH, MESH)

color_axis()
#
#grid(P(-MAX,-MAX,-MAX), P( MAX, MAX,-MAX), 4, 4)
#h_axis_labels(P(-MAX, MAX,-MAX), P(MAX, MAX,-MAX), 2*MAX, P( 2,-2), br)
#

label(P(0,0.25+MAX,0), P(4,0), "$y$", r)
label(P(0,0,0.25+MAX), P(0,4), "$z$", t)

# back 1/2 of surface (x<0), drawn back-to-front
clip()
fill()
for i in range(-N, 0):
  for j in range(-N, N):
    draw_mesh_quad(i,j)

fill(false)
clip(false)

# coordinate axes
bold()
use_pstricks()
ps_axis_color()
psset("linecolor=axis_color,fillcolor=axis_color")
#  color_axis();
dart(P(0,-MAX,0), P(0,0.25+MAX,0))
dart(P(0,0,0), P(0,0,0.25+MAX))
use_pstricks(false)

# vertices of slicing planes
pt1 = P(x_0, 0, z_0)
pt2 = P(x_0, 0, -MAX)
pt3 = P(x_0, MAX, -MAX)
pt4 = P(0, y_0, z_0)
pt5 = P(0, y_0, -MAX)
pt6 = P(MAX, y_0, -MAX)
pt7 = P(x_0, 0, 0)
pt8 = P(x_0, MAX, z_0)
pt9 = P(0, y_0, f(0,y_0))
pt10 = P(MAX, y_0, z_0)

# lower outlines of slicing planes (partially hidden by surface)
color_xslice()
polyline((pt1, pt2, pt3)).draw()
color_yslice()
polyline((pt4, pt5, pt6)).draw()

plain()
clip()
fill()
# front 1/2 of surface (x>0)
for i in range(0,N):
  for j in range(-N,N):
    draw_mesh_quad(i,j)

fill(false)
clip(false)

# last coordinate axis
bold()
use_pstricks()
psset("linecolor=axis_color, fillcolor=axis_color")
dart(P(-MAX,0,0), P(0.25+MAX,0,0))
use_pstricks(false)

color_axis()
label(P(0.25+MAX,0,0), P(-2,-2), "$x$", bl)

# slices of the graph and labels
pen(1.5)
color_yslice()
plot_triple(Df_x, -x_0, MAX-x_0, 20)
label(Df_x(MAX-x_0), P(-4,0), 
"\\scriptsize $\\displaystyle\\frac{\\partial f}{\\partial x}$: $y$ constant", l)

clip()
color_xslice()
plot_triple(Df_y, -y_0, MAX-y_0, 20)
clip(false)
label(Df_y(MAX-y_0), P(4,0), 
"\\scriptsize $\\displaystyle\\frac{\\partial f}{\\partial y}$: $x$ constant", br)

# overall label
black()
label(P(-MAX, 0, MAX), P(2,2), 
"\\scriptsize $z=\\displaystyle\\frac{1}{2}(y^3-3x^2y)$", tr)

# fake transparency: redraw slicing planes thinly
clip()
pen(0.15)
color_xslice()
polyline((pt1, pt2, pt3)).draw()
color_yslice()
polyline((pt4, pt5, pt6)).draw()

# upper outlines of slicing planes (not hidden by surface) 
bold()
color_xslice()
polyline((pt7, pt1, pt8, pt3)).draw()

color_yslice()
polyline((pt9, pt4, pt10, pt6)).draw()

end()