## -*-SM-*- Set SM-mode in emacs
# This file contains all SM macros for the user interface of SURFACE:
#    surface viewpoint box3 label3 projection axis3 getticks trunc minmaxlim
#
# To overload SURFACE and VIEWPOINT say
#	overload viewpoint 1
#	overload surface 1
#
# To try a sample image, use
#
#   testsurface
#
# To see all the SM variables related to Surface, Viewpoint and box3, use
#
#   Show3
#
# Surface needs no parameters, and if it has it is very flexible about them.
# Viewpoint does the same as VIEWPOINT but stores the parameters.
# box3 draws 3 projected axes. In general use it after calling surface.
# label3 puts a label to the axis given in its first parameter (=x,y, or z).
# projection can be used to project 3d coordinates to 2d coordinates.
# axis3 draws a projected axis with endpoints given in 3d coordinates.
# getticks can be used to get good ticksize values to avoid crowding.
# trunc is the mathematical int function.
# minmaxlim gets min & max of the image within given limits and restores limits
testsurface	# This is an example for using Surface
		image(20,20)
		set kx=0.1*pi set ky=0.2*pi
		set ixs=0,19
		do iy=0,19{
		   set image[ixs,$iy]=cos(ixs*kx+1)*sin($iy*ky+1)
		}
		if(!$?Theta) { define Theta 10 }
		if(!$?Phi) { define Phi 20 }
		if(!$?Distance) { define Distance 0 }
		define Theta ?
		define Phi ?
		define Distance ?
		Viewpoint $Theta $Phi $Distance
		Surface
		box3
		label3 x X-label [xunit]
		label3 y Y-label [yunit]
		label3 z Z-label [zunit]
Show3		# Show all variables related to surface drawing. We echo
		# the actual variable names so the user may use them.
		echo Theta_view=$Theta_view Phi_view=$Phi_view\
		    Distance_view=$Distance_view  Obey_limits=$?Obey_limits
		echo Lo_x=$Lo_x Hi_x=$Hi_x Lo_y=$Lo_y Hi_y=$Hi_y Lo_z=$Lo_z\
		    Hi_z=$Hi_z
		foreach 0 {x y z} {
		   if($?Label_x_$0) {
		      echo Label_x_$0=$Label_x_$0 Label_y_$0=$Label_y_$0\
		 Label_a_$0=$Label_a_$0 Label_p_$0=$Label_p_$0
		   }
		}
surface		# Overloaded version of SURFACE
		Surface
Surface	05	# Usage: Surface [type] [z1 z2] [x y]
		# where x and y must be arrays of dimen > 1, and
		# type, z1 and z2 are numbers.
		# If type is given it has to be the first parameter.
		# Passes parameters (defaults if necessary) for SURFACE and 
		# defines Lo_x Hi_x Lo_y Hi_y Lo_z Hi_z Obey_limits for box3.
		define NX image
		define NY image
		if(!$?NX) {
		   echo No image is defined
		   return
		}
		define _Nparam 0
		do _P=1,5 {
		   if($?$_P) { define _Nparam $_P }
		}
		if($_Nparam==1 || $_Nparam==3 || $_Nparam==5) {
		   define Type_surface $1
		   define _Pmin 2
		} else {
		   define Type_surface 3
		   define _Pmin 1
		}
		if(int($Type_surface/10)-10*int($Type_surface/100)) {
		   define Obey_limits 1
		} else {
		   define Obey_limits DELETE
		}
		do _P=$_Pmin,$_Nparam {
		   set _tmp=$$_P
		   if (dimen(_tmp)>1) {
		      if(!$?X_done) {
		         set x_surface=$$_P define X_done yes
		      } else {
		         set y_surface=$$_P define Y_done yes
		      }
		   } else {
		      if(!$?Lo_done) {
		         define Lo_z $$_P define Lo_done yes
		      } else {
		         define Hi_z $$_P define Hi_done yes
		      }
		   }
		}
		if($?Y_done) {
	#Note: due to a bug MINMAX doesn't work if $fx1>$fx2 or $fy1>$fy2
	#but here we make sure that Lo_x<Hi_x and Lo_y<Hi_y
		   define Lo_x (x_surface[0])
		   define Hi_x (x_surface[dimen(x_surface)-1])
		   if($Lo_x > $Hi_x) {
		      define _Tmp $Lo_x define Lo_x $Hi_x define Hi_x $_Tmp
		   }
		   define Lo_y (y_surface[0])
		   define Hi_y (y_surface[dimen(y_surface)-1])
		   if($Lo_y > $Hi_y) {
		      define _Tmp $Lo_y define Lo_y $Hi_y define Hi_y $_Tmp
		   }
		} else {
		   define Lo_x 0
		   define Hi_x ($NX-1)
		   define Lo_y 0
		   define Hi_y ($NY-1)
		}
		if (!$?Hi_done) {
		   if($Type_surface<100) {
		      minmaxlim $Lo_x $Hi_x $Lo_y $Hi_y Lo_z Hi_z
		   } else {
		      minmaxlim 0 $NX 0 $NY Lo_z Hi_z
		   }
		}
		if($Lo_z==$Hi_z) {
		   set _tmp=($Lo_z)? 0.05*abs($Lo_z):0.05
		   define Lo_z ($Lo_z-_tmp) define Hi_z ($Hi_z+_tmp)
		}
		if ($?Y_done) {
		   SURFACE $Type_surface $Lo_z $Hi_z x_surface y_surface
		} else {
		   SURFACE $Type_surface $Lo_z $Hi_z
		}
		foreach 0 {Type_surface NX NY _Pmin _Nparam Lo_done Hi_done \
		   X_done Y_done _Tmp} {define $0 DELETE}
		foreach 0 {x_surface y_surface _tmp}{DELETE $0}
box3		# Usage: box3
		# Draw a 3D-box defined by Lo_x Hi_x Lo_y Hi_y Lo_z Hi_z and
		# Phi_view Theta_view Distance_view, and according to
		# Obey_limits. These are set by Viewpoint and Surface, but 
		# box3 has defaults for them. Axis-label variables are set. 
		if(!$?Lo_x || !$?Hi_x || !$?Lo_y || !$?Hi_y || \
		    !$?Lo_z ||!$?Hi_z) {
		   define NX image define NY image
		   if(!$?NX) {
		      define Lo_x 0 define Hi_x 1
		      define Lo_y 0 define Hi_y 1
		      define Lo_z 0 define Hi_z 1
		   } else {
		      define Lo_x 0 define Hi_x ($NX-1)
		      define Lo_y 0 define Hi_y ($NY-1)
		      minmaxlim 0 $NX 0 $NY Lo_z Hi_z
		      if($Lo_z==$Hi_z) {
		         set _tmp=($Lo_z)? 0.05*abs($Lo_z):0.05
		         define Lo_z ($Lo_z-_tmp) define Hi_z ($Hi_z+_tmp)
		      }
		   }
		}
		if(!$?Phi_view || !$?Theta_view || !$?Distance_view) {
		   define Phi_view (90/pi)
		   define Theta_view (90/pi)
		   define Distance_view 0
		}
		define Fr_y ((abs($Phi_view)<90)? $Lo_y : $Hi_y)
		define Bk_y  ((abs($Phi_view)<90)? $Hi_y : $Lo_y)
		set edge_x=<$Lo_x $Hi_x $Lo_x $Hi_x $Lo_x $Hi_x $Lo_x $Hi_x>
		set edge_z=<$Lo_z $Lo_z $Hi_z $Hi_z $Lo_z $Lo_z $Hi_z $Hi_z>
		set edge_y=<$Fr_y $Fr_y $Fr_y $Fr_y $Bk_y $Bk_y $Bk_y $Bk_y>
		projection edge_x edge_y edge_z edge_px edge_py
		set back_axial= $Distance_view<0 && abs($Phi_view)>90
		if(!$?Obey_limits) {    #find limits for the projected 3D box 
		                        # from edge_px and edge_py
		   foreach _X {px py} {
		      set edge_=edge_$_X
		      sort <edge_>
		      define Lo_$_X (edge_[0])
		      define Hi_$_X (edge_[7])
		      set _tmp=($Lo_$_X<$Hi_$_X)? $Hi_$_X-$Lo_$_X : \
		      		      ( ($Lo_$_X)? abs($Lo_$_X) : 1 )
		      define Lo_$_X ($Lo_$_X-0.05*_tmp)
		      define Hi_$_X ($Hi_$_X+0.05*_tmp)
		   }
		   if(back_axial) {
		      define _Tmp $Lo_px define Lo_px $Hi_px define Hi_px $_Tmp
		   }
		   if( $Lo_z > $Hi_z ) {
		      define _Tmp $Lo_py define Lo_py $Hi_py define Hi_py $_Tmp
		   }
		   limits $Lo_px $Hi_px $Lo_py $Hi_py
		}
		#Draw the 3 axes
		if($Theta_view>=0) {
		   define _Z $Lo_z define _C 2
		} else {
		   define _Z $Hi_z define _C 3
		}
		if(abs($Phi_view)<90) {
		   axis3 x $Lo_x $Fr_y $_Z $Hi_x $Fr_y $_Z 1 $_C
		} else {
		   axis3 x $Hi_x $Fr_y $_Z $Lo_x $Fr_y $_Z 1 $_C
		}
		if((edge_px[0]<=edge_px[4])!=back_axial) {
		   define _X ((abs($Phi_view)<90)? $Hi_x : $Lo_x)
		   axis3 y $_X $Fr_y $_Z $_X $Bk_y $_Z 1 $_C
		} else {
		   define _X ((abs($Phi_view)>90)? $Hi_x : $Lo_x)
		   axis3 y $_X $Bk_y $_Z $_X $Fr_y $_Z 1 $_C
		}
		define _X (((edge_px[0]< edge_px[1])==back_axial)? $Hi_x: $Lo_x)
		define _Y (((edge_px[0]<=edge_px[4])!=back_axial)? $Fr_y: $Bk_y)
		axis3 z $_X $_Y $Lo_z $_X $_Y $Hi_z 2 5
		foreach 0 < edge_x edge_y edge_z edge_px edge_py edge_ _tmp \
		    back_axial > {
		   DELETE $0}
		foreach 0 { Fr_y Bk_y Lo_px Lo_py Hi_px Hi_py _Tmp _X _Y _Z _C \
		    NX NY } {
		   define $0 DELETE}
viewpoint	# Overloaded version of VIEWPOINT
		Viewpoint
Viewpoint 3	# stores parameters of VIEWPOINT for box3, check for bad angles
		if($1<-89||$1>89) {
		   echo Please choose theta in [-89,89]
		   return
		}
		if($2<-179.5||$2>179.5) {
		   echo Please choose phi in [-179.5,179.5]
		   return
		}
		define Theta_view $1
		define Phi_view   $2
		define Distance_view $3
		if($2 >-90.5 && $2 <=-90) {define Phi_view -90.5}
		if($2 <-89.5 && $2 > -90) {define Phi_view -89.5}
		if($2 > -0.5 && $2 <=  0) {define Phi_view  -0.5}
		if($2 <  0.5 && $2 >   0) {define Phi_view   0.5}
		if($2 > 89.5 && $2 <= 90) {define Phi_view  89.5}
		if($2 < 90.5 && $2 >  90) {define Phi_view  90.5}
		VIEWPOINT $Theta_view $Phi_view $Distance_view
projection 5	# Usage: projection x,y,z,xprojected,yprojected
		# uses Lo_x Hi_x Lo_y Hi_y Lo_z Hi_z Phi_view Theta_view
		# Distance_view (0 for infinite projective, and negative
		# for axonometric projection)
		set front_y=(abs($Phi_view)<90)? $Lo_y : $Hi_y
		set _phi=$Phi_view/180*pi set _theta=$Theta_view/180*pi
		if($Hi_z < $Lo_z) { 
		   set _theta = -_theta 
		}
		if($Distance_view < 0) {
		   set _diy= (abs(_phi)<=pi/2)? 1 : -1
		   set $4=$1-_diy*tan(_phi)*\
		       abs(($2-front_y)*($Hi_x-$Lo_x)/($Hi_y-$Lo_y))
		   set $5=$3+tan(_theta)*\
		       abs(($2-front_y)*($Hi_z-$Lo_z)/($Hi_y-$Lo_y))
		   foreach 0 {_diy _phi _theta front_y}{DELETE $0}
		   return
		}
		set front_x=($Phi_view>0)? $Lo_x : $Hi_x
		set front_z=($Theta_view<0)? $Lo_z : $Hi_z
		set x_=($1-front_x)/abs($Hi_x-$Lo_x)
		set y_=($2-front_y)/abs($Hi_y-$Lo_y)
		set z_=($3-front_z)/abs($Hi_z-$Lo_z)
		set $4=cos(_phi)*x_ - sin(_phi)*y_
		set $5=sin(_phi)*sin(_theta)*x_ + \
		    cos(_phi)*sin(_theta)*y_+cos(_theta)*z_
		if($Distance_view>0) {
		   set _s_=1+(sin(_phi)*cos(_theta)*x_ + \
		       cos(_phi)*cos(_theta)*y_ -sin(_theta)*z_)/$Distance_view
		   set $4=$4/_s_
		   set $5=$5/_s_
		}
		foreach 0 {front_x front_y front_z x_ y_ z_ _phi _theta _s_} {
		   DELETE $0
		}
axis3 9		# Usage: axis3 x/y/z,x_1,y_1,z_1,x_2,y_2,z_2,ilabel,iclock
		# The first parameter is the name of the axis (x, y or z),
		# the endpoints are in 3D coordinates, and the axis is 
		# projected to 2D. The number of big ticks is set to 2, 
		# number of smallticks are determined by getticks.
		# Sets axis-label screen coordinates: Label_x_$1, Label_y_$1;
		# the angle: Label_a_$1; and the putlabel parameter: Label_p_$1.
		set x_1=$2 set y_1=$3 set z_1=$4
		set x_2=$5 set y_2=$6 set z_2=$7
		set a_1=$1_1 set a_2=$1_2
		set big_=2
	# Due to some bug, axis does not work with _lab if a_1>a_2, we fix it
		if(a_1>a_2) {set a_1=-a_1 set a_2=-a_2}
		getticks a_1 a_2 sml_ big_
		set tmp_=trunc($(a_1/sml_)),(trunc($(a_2/sml_)))
		set sml_=tmp_*sml_
		set tmp_=trunc($(a_1/big_)),(trunc($(a_2/big_)))
		set big_=tmp_*big_
		set lab_=($1_1>$1_2)? sprintf('%g',-big_) : sprintf('%g',big_)
		projection x_1 y_1 z_1 px_1 py_1
		projection x_2 y_2 z_2 px_2 py_2
		set a_21=a_2-a_1
		set tmp_x=sml_/a_21*(x_2-x_1)+a_2/a_21*x_1-a_1/a_21*x_2
		set tmp_y=sml_/a_21*(y_2-y_1)+a_2/a_21*y_1-a_1/a_21*y_2
		set tmp_z=sml_/a_21*(z_2-z_1)+a_2/a_21*z_1-a_1/a_21*z_2
		projection tmp_x tmp_y tmp_z sml_px sml_py
		set tmp_x=big_/a_21*(x_2-x_1)+a_2/a_21*x_1-a_1/a_21*x_2
		set tmp_y=big_/a_21*(y_2-y_1)+a_2/a_21*y_1-a_1/a_21*y_2
		set tmp_z=big_/a_21*(z_2-z_1)+a_2/a_21*z_1-a_1/a_21*z_2
		projection tmp_x tmp_y tmp_z big_px big_py
		if(px_2==px_1) {
		   set py_21=py_2-py_1
		   set sml_=(py_2-sml_py)/py_21*a_1+(sml_py-py_1)/py_21*a_2
		   set big_=(py_2-big_py)/py_21*a_1+(big_py-py_1)/py_21*a_2
		} else {
		   set px_21=px_2-px_1
		   set sml_=(px_2-sml_px)/px_21*a_1+(sml_px-px_1)/px_21*a_2
		   set big_=(px_2-big_px)/px_21*a_1+(big_px-px_1)/px_21*a_2
		}
		set gx_1=int((px_1*($gx2-$gx1)+$gx1*$fx2-$gx2*$fx1)/($fx2-$fx1))
		set gx_2=int((px_2*($gx2-$gx1)+$gx1*$fx2-$gx2*$fx1)/($fx2-$fx1))
		set gy_1=int((py_1*($gy2-$gy1)+$gy1*$fy2-$gy2*$fy1)/($fy2-$fy1))
		set gy_2=int((py_2*($gy2-$gy1)+$gy1*$fy2-$gy2*$fy1)/($fy2-$fy1))
		set g_len=int(sqrt((gx_1-gx_2)**2+(gy_1-gy_2)**2))
		angle $( 180/pi*atan2($((gy_2-gy_1)*$aspect),$(gx_2-gx_1)) )
		axis $(a_1) $(a_2) sml_ big_ lab_ "$!(gx_1)" "$!(gy_1)" \
		    $(g_len) $8 $9
		#Find distance for axis-label from the midpoint
		if($8==1) {
		   # For labels parallel to the axis, use string height
		   set tmp_=length(lab_[0])
		   define sheight |
		   set label_dist=$sheight
		} else {
		   # For labels orthogonal to the axis, find the longest label
		   set tmp_=-length(lab_)
		   sort{tmp_}
		   set label_dist=-tmp_[0]*$aspect
		}
		# Depending on the parity of ICLOCK the axis-label is put above
		# or below the axis. Aspect ratio is taken into account.
		set angle_ort=pi/180*( (int($9/2)==$9/2)? $angle-90: $angle+90 )
		set cos_o=cos(angle_ort)*$aspect
		set sin_o=sin(angle_ort)
		# LABEL_OFFSET=400 from axis.c, thus 2.5*400=1000 will do
		set label_dist=(label_dist+$expand*1000)/sqrt(cos_o**2+sin_o**2)
		define Label_x_$1 $((gx_1+gx_2)/2+label_dist*cos_o)
		define Label_y_$1 $((gy_1+gy_2)/2+label_dist*sin_o)
		define Label_a_$1 $angle
		# The putlabel parameter is 2 for anticlockwise, 8 for clockwise
		define Label_p_$1 $((int($9/2)==$9/2)? 2: 8)
		angle 0
		foreach 0 {x_1 x_2 px_1 px_2 px_12 y_1 y_2 py_1 py_2 py_12 \
		    z_1 z_2 a_1 a_2 a_21 big_ big_px big_py sml_ sml_px sml_py \
		    tmp_ tmp_x tmp_y tmp_z gx_1 gx_2 gy_1 gy_2 g_len \
		    label_dist angle_ort cos_o sin_o} {
		   DELETE $0
		}
	#Due to another mysterious bug, lab_ should be deleted separately,
	#otherwise the HISTORY file may be deleted (e.g. when lab_ contains "0")
		DELETE lab_
		define sheight DELETE
label3 122	#Usage: label3 x/y/z label_text
		if('$1'!='x'&&'$1'!='y'&&'$1'!='z') {
		   echo First parameter should be x, y, or z
		   return
		}
		if(!$?Label_x_$1) {
		   echo Label position is not defined, use box3 to set it
		   return
		}
		define 0 $angle
		angle $Label_a_$1
		RELOCATE ($Label_x_$1 $Label_y_$1)
		putlabel $Label_p_$1 $2
		#putlabel 5 \raise$Label_d_$1"$!2"
		angle $0
getticks 4	# xmin xmax 0 number_of_big_ticks --> xmin xmax smallx bigx
		set _smallticks={1 0.1 0.2 0.5 0.5 0.5 0.5 0.5 1 0.5}  
		set _range=abs($2 - $1)
		set _nbig=$4
		set _order=10**trunc($(lg(_range)))
		set _norder=trunc($(_range/_order))
		set $3=_order*_smallticks[_norder]
		set $4=_order*_norder/_nbig
		foreach 0 {_smallticks _range _nbig _order _norder}{
		   DELETE $0
		}
trunc	1	# Calculates mathematical int function of $1
		set $0 = ($1 >= int($1)) ? int($1) : int($1) - 1
minmaxlim 6	# xmin,xmax,ymin,ymax,imagemin,imagemax
		# same as limits and minmax but original limits are restored
		foreach 0 {fx1 fx2 fy1 fy2} {
		   define old_$0 $$0
		}
		limits $1 $2 $3 $4
		minmax $5 $6
		limits $old_fx1 $old_fx2 $old_fy1 $old_fy2
		foreach 0 {fx1 fx2 fy1 fy2} {
		   define old_$0 DELETE
		}