*! version 1.2.0 , apr 98, Guy van Melle  STB-45 gr30
*!
*! 3D plot with hidden-line-removal
*!
*! syntax:	hidlin  func  
*!		[,  x(lo hi step)  y(lo hi step)  Eye(xE yE zE)
*!		Neg  Box Coo Lines([x][y])  Text  TMag(int 100) 
*!		XMar(int 10)  YMar(int 10)]  SAVing(filename) ]
*!
*! function func z=f(x,y) must be declared in hlfunc.ado
*
*--------------
prog def hidlin
*--------------
	version 5.0
	if "`*'"=="clean" { * erase macs and mats
		hlclean
		exit
	}		
	qui des,sh
	if _result(1) {
		di in blu "warning: data in memory will be lost"
		di in blu "OK to continue (Y/N) ? " _req(hlYesNo)
		if substr(upper("$hlYesNo"),1,1)!="Y" { exit }
		glo hlYesNo
	}
	drop _all

*  keep it clean if break; don't allow more than 500*500 grid

	preserve   
	loc a=real("0$hlNmax")
	if `a'==. | `a'==0 | `a'>1000 { glo hlNmax 1001 }
	
	if "`*'"!="" {
		parse "`*'", parse(",") 
		if "`1'"!="," {
			loc p `1'
			glo hlprog "hl`p'"
			confirm file "${hlprog}.ado"
			mac shift
		}
		Get_info `*'
		if "$S_1"!="" {	
			confirm new file $S_1
			loc sav , saving($S_1) 
		}
		glo S_1
	}
	projmat
	loc xmx= 32000 - 3
	loc ymx= 23060 - 1500 * (($hlBox | $hlCoo) + $hlTex)
	winsiz `xmx' `ymx'
	base
	loc p=cond($hlPen,$hlPen,1)

	gph open `sav'
		gph pen `p'
		gph vline gy gx			
		if $hlXli {  xlin }
		if $hlYli {  ylin }
		finish `xmx' `ymx' 
		if $hlTex { hltex }
	gph  close
end
*
*----------------
prog def Get_info
*
	#delimit ;
	loc options	"Xinfo(string) Yinfo(string) Eye(string) 
			Lines(string) Box Coo Neg Penc XMar(int 10) 
			YMar(int 10) Text TMag(int 100) SAVing(string)" 
	; #delimit cr

	parse "`*'"
	if "`xinfo'"=="" { loc xinfo $hlX }
	else {
		info 1 `xinfo'
		loc xinfo $S_1
		glo  hlX  $S_1
	}
	if "`yinfo'"=="" { loc yinfo $hlY }
	else {
		info 1 `yinfo'
		loc yinfo $S_1
		glo  hlY  $S_1
	}
	if  "`eye'" =="" { loc  eye  $hlE }
	else {
		info 0  `eye'
		loc  eye  $S_1
		glo  hlE  $S_1
	}
	if "`xinfo'"=="" | "`yinfo'"=="" | "`eye'"=="" {
		di in red "xinfo, yinfo, eye must  all  be provided"
		di in blu "x(lo hi step)  y(lo hi step)  eye(x y z)"
		exit 198
	}
	if "`lines'"!="" { 
		glo hlXli=index("`lines'","x")>0 
		glo hlYli=index("`lines'","y")>0 
	} 
	else {
		glo hlXli= real("0$hlXli")
		if $hlXli== .  { glo hlXli=0 }
		glo hlYli= real("0$hlYli")
		if $hlYli== .  { glo hlYli=0 }
	}
	if $hlXli+$hlYli==0 {
		glo hlXli=1
		glo hlYli=1
	}
	glo hlBox=     ("`box'" !="")
	glo hlCoo=     ("`coo'" !="")
	glo hlNeg=     ("`neg'" !="")
	glo hlTex=     ("`text'"!="")
	glo hlMag= 	 `tmag'
	glo hlYma= 	 `ymar'
	glo hlXma= 	 `xmar'
	glo hlPen= cond("`penc'"=="",2,0)
	glo S_1		 `saving'
end
*------------
prog def info
*
	loc veri `1'
	mac shift
	loc info `*'
	parse   "`*'", parse(" ,")
	loc k 1 
	while "``k''"!="" { 
		if "``k''"=="," { loc `k' " "}
		else            { loc `k'=``k''} /*allows for expr in info*/
	loc  k=`k'+1
	}
	loc inf `*'
	loc w: word count `inf'
	loc ok=`w'==3
	if `ok' & `veri' {
		loc a: word 1 of `inf'
		loc b: word 2 of `inf'
		if `a'>`b' {
			loc c `a'
			loc a `b'
			loc b `c'
		}
		loc c: word 3 of `inf'
		loc c= abs(`c')
		if `c'<`b'-`a'	{ loc inf "`a' `b' `c'" }
		else			{ loc ok 0 }
	}
	if `ok' { glo S_1 `inf' }
	else { 
		di in red "`info' invalid"
		exit 198
	}
end
*---------------
prog def projmat
*
	tempname e1 e2 e3 d2 d3 r

	loc e1: word 1 of $hlE
	loc e2: word 2 of $hlE
	loc e3: word 3 of $hlE

	if `e1'*`e2'==0 {
		di in red "you can't sit there"
		exit 198
	}
	sca `d2'= sqrt(`e1'*`e1' + `e2'*`e2')
	sca `d3'= sqrt(`d2'*`d2'+`e3'*`e3')

	sca `r' =   `e3'/`d3'
	loc  a  = - `e2'/`d2'
	loc  b  =   `e1'/`d2'
	loc  d  = - `b' * `r'
	loc  e  =   `a' * `r'
	loc  f  =   `d2'/`d3' 

	mat  hlPj= (`a' `d' \ `b' `e' \ 0 `f')
end
*--------------
prog def projec
*
* call: projec [in L1/L2]  [, Drop] 
*
	loc in "opt"
	loc options "Drop"
	parse "`*'" 
	tempvar  z

	$hlprog x y `z' `in'
	qui {
		if $hlNeg { replace `z'= - `z' }	
		if "`drop'"!="" {
			cap drop px
			cap drop py
			g px=.
			g py=.
		}
		replace px= x*hlPj[1,1] + y*hlPj[2,1] + `z'*hlPj[3,1] `in'
		replace py= x*hlPj[1,2] + y*hlPj[2,2] + `z'*hlPj[3,2] `in'
		if "`drop'"=="" {
			replace gx= round(px * $hlXfa + $hlXsh, 1) `in'
			replace gy= round(py * $hlYfa + $hlYsh, 1) `in'
		}
 	}
end
*--------------
prog def winsiz
*
	loc xmx= `1'
	loc ymx= `2'

	loc a "( $hlX \ $hlY )"
	mat hlWk= `a'
	loc r1 = hlWk[1,2]-hlWk[1,1]
	loc r2 = hlWk[2,2]-hlWk[2,1]

	loc d1= abs(hlWk[1,3])
	loc d2= abs(hlWk[2,3])

	loc e1: word 1 of $hlE
	loc e2: word 2 of $hlE

	loc e1= abs(`e1')
	loc e2= abs(`e2')

	loc c = (`r1'/`e1'+`r2'/`e2')/(`r1'/`d1'+`r2'/`d2')

	loc f1=`e1'*`c'
	loc f2=`e2'*`c'

	loc ok 0
	while !`ok' {
		glo hlXn= int(abs(`r1'/`f1'+0.01))
		glo hlYn= int(abs(`r2'/`f2'+0.01))
		glo hln = $hlXn + $hlYn + 1
        loc  ok = $hln <= $hlNmax
		if !`ok' {
			loc f1=`f1'* $hln / $hlNmax
			loc f2=`f2'* $hln / $hlNmax
		}
	}
	mat hlWk[1,3]=`f1'
	mat hlWk[1,2]= hlWk[1,1]+$hlXn*`f1'
	mat hlWk[2,3]=`f2'
	mat hlWk[2,2]= hlWk[2,1]+$hlYn*`f2'

	tempname U crn
	mat `U'=hlWk
	sca `crn'=0
	while `crn'!= 1 { 
		orient `crn' 
		if `crn'>2 {           
			mat  `U'[1,1]=  hlWk[1,2]
			mat  `U'[1,2]=  hlWk[1,1]
			mat  `U'[1,3]= -hlWk[1,3]
		}
		if mod(`crn',2)==0 { 
			mat  `U'[2,1]= hlWk[2,2]
			mat  `U'[2,2]= hlWk[2,1]
			mat  `U'[2,3]=-hlWk[2,3]
		}
		if `crn'!=1 { mat hlWk=`U' }
	}
	mat hlCx= J(1,4,0)
	mat hlCy= J(1,4,0)
	mat hlCz= J(1,4,0)
	loc c 0
	while `c'<4 {
		loc c=`c'+1
		mat hlCx[1,`c']= px[`c']
		mat hlCy[1,`c']= py[`c']
		mat hlCz[1,`c']=  z[`c']
	}
* this is the slow part so reduce steps temporarily if necessary
	
	loc div= round($hln/200+.499,1)
	loc nx = int( ($hlXn+1)/`div')
	loc ny = int( ($hlYn+1)/`div')

	drop _all
	qui {	
		set obs `nx'
		g x =_n-1
		expand  `ny'
		sort x
		g y = mod(_n-1,`ny')

		replace x = hlWk[1,1] + x * hlWk[1,3] * `div'
		replace y = hlWk[2,1] + y * hlWk[2,3] * `div'
	} 
 	projec , d 

	qui su px
	loc a1= _result(5)
	loc a2= _result(6)
	loc a = `a2'-`a1'
	loc ma= $hlXma/100
   
	glo hlXfa =        (`xmx'*( 1 - 2 * `ma') / `a')
	glo hlXsh =   round(`xmx'*( `ma' *(`a2'+`a1') - `a1')/`a', 1)

	qui su py
	loc b1= _result(5)
	loc b2= _result(6)
	loc b = `b2'-`b1'
	loc mb= $hlYma/100

	glo hlYfa = -      (`ymx'*( 1 - 2 * `mb') / `b')
	glo hlYsh = - round(`ymx'*( `mb' *(`b1'+`b2') - `b2')/`b', 1)

	loc j=0
	while `j'<4 {
		loc j=`j'+1
		mat hlCx[1,`j']= round( hlCx[1,`j'] * $hlXfa + $hlXsh, 1)
		mat hlCy[1,`j']= round( hlCy[1,`j'] * $hlYfa + $hlYsh, 1)
	}
end
*--------------
prog def orient
*
	qui {  
		drop _all
		set obs 4  /*(hlWkx1,y1) (Wkx1,y2) (Wkx2,y1) (Wkx2,y2)*/
		g x=.
		g y=.
		g z=.
		loc c 0
		loc i 0
		while `i'<2 {
			loc i=`i'+1
			loc j 0
			while `j'<2 {
				loc j=`j'+1
				loc c=`c'+1
				replace x = hlWk[1,`i'] in `c'
				replace y = hlWk[2,`j'] in `c'
			}
		}
	}
	projec  , d

	qui replace z=(py-x*hlPj[1,2]-y*hlPj[2,2])/hlPj[3,2]
	loc a 1
	loc c 1
	loc b py[1]
	while `a'<4 {
		loc a =`a'+1
		loc p = py[`a']
		if `p'<`b' {
			loc c=`a'
			loc b=`p'
		}
	}
	sca `1'=`c' /* lowest proj.corner */
end  
*------------
prog def base
*
	drop _all

	loc ny= $hlYn+1
	qui {
		set obs $hln

		g int 	n =_n-`ny'
		g		x = hlWk[1,1]					in  1/`ny'
		replace x = x[_n-1] + hlWk[1,3]			if _n>`ny'

		g		y = hlWk[2,2]-(_n-1)*hlWk[2,3]	in  1/`ny'
		replace y = hlWk[2,1]					if _n>`ny'

		projec , d

		g int  gx = round(px * $hlXfa + $hlXsh, 1)
		g int  gy = round(py * $hlYfa + $hlYsh, 1) 
		g int  ba= gy
		g int  lo= .
		g int  hi= .
		g byte isVi=3
	}
end
*--------------
prog def finish
*
	loc xmx `1'
	loc ymx `2'
  
	if $hlBox | $hlCoo {
		if (hlCx[1,1]-hlCx[1,2]) * (hlCx[1,1]-hlCx[1,3]) < 0 {
			loc a 1,1
			loc b 1,2
			loc c 1,3
			loc d 1,4
		}
		else {
			loc q= 2 + (hlCy[1,2] < hlCy[1,3])
			loc a 1,`q'
			loc b 1,1
			loc c 1,4
			loc q=5-`q' 
			loc d 1,`q'
		}
		loc y1= hlCy[`a']
		loc y2= hlCy[`b']
		loc y3= hlCy[`c']
		loc y4= hlCy[`d']
		loc x1= hlCx[`a']
		loc x2= hlCx[`b']
		loc x3= hlCx[`c']
		loc x4= hlCx[`d']
   	}
*--- show box ---

	if $hlBox {	
		loc dv= abs( $hlYfa * hlPj[3,2] ) 
		loc dd= max(`y1'+hlCz[`a']*`dv',`y2'+hlCz[`b']*`dv',`y3'+hlCz[`c']*`dv')
*		loc dd= max(hlCz[`a']*`dv',hlCz[`b']*`dv',hlCz[`c']*`dv')*.9

		loc y5= `y1'+.85*(`ymx'+hlCz[`a']*`dv'-`dd')
		loc y6= `y2'+.85*(`ymx'+hlCz[`b']*`dv'-`dd')
		loc y7= `y3'+.85*(`ymx'+hlCz[`c']*`dv'-`dd')

		gph line `y1' `x1' `y5' `x1'
		gph line `y2' `x2' `y6' `x2'
		gph line `y3' `x3' `y7' `x3'
		gph line `y5' `x1' `y6' `x2'  
		gph line `y5' `x1' `y7' `x3'  
    }
*--- show corner coo ---
 
	if $hlCoo{	
		if $hlPen { gph pen 1}
		gph font 520 260
		loc al = 300
		loc du = 800
		loc j 0
		while `j'<4 {
			loc j = `j'+1
			loc q = substr("abcd",`j',1)
			loc t = substr("``q''", 3 ,1)
			loc dv= 2*( 2*`x`j'' < `xmx')-1

			loc u = `y`j''- 3*`du'*(`j'==4)
			loc v = `x`j''-`al'*`dv'
                
			loc u = `u'+`du'                
			loc r : display %6.2f scalar(hlWk[1,1+(`t'>2)])
			gph text `u' `v' 0  `dv' `r'

			loc u = `u'+`du'
			loc r : display %6.2f scalar(hlWk[2,1+(mod(`t',2)==0)])
			gph text `u' `v' 0  `dv' `r'

			loc u = `u'+`du'
			loc r : display %6.2f scalar(hlCz[``q''])
			gph text `u' `v' 0  `dv' `r'
		}
		if $hlPen { gph pen $hlPen}
	}
	else {
*		sho arrows: from lox dir hix & from loy dir hiy
*		not implemented
	}
end
*------------
prog def xlin
*
	qui {
		replace lo= ba
		replace hi= lo
		loc L1= $hlYn+1
		loc L2= $hln
		replace x= hlWk[1,1] + (_n-`L1') * hlWk[1,3] in `L1'/`L2'
	}
	loc k  0
	while `k'<= $hlYn {
		qui replace y= hlWk[2,1] +   `k' * hlWk[2,3] in `L1'/`L2'

		projec in `L1'/`L2' 
        	show	  `L1' `L2'

		loc L1=`L1'-1
		loc L2=`L2'-1
		if `L1'> 0  { qui replace x = x[_n+1]  in `L1'/`L2'}
		loc k =`k' +1
	} 
end
*------------
prog def ylin
*
/*
  neat: y-lines are drawn the other way around, so switch 
  base-curve upside down and use exactly same logic as for x-lines
*/
	qui {
		replace lo= ba[_N+1-_n] 
		replace hi= lo
		loc L1= $hlXn+1
		loc L2= $hln
		replace y= hlWk[2,1] + (_n-`L1') * hlWk[2,3] in `L1'/`L2'
	}
	loc k  0
	while `k'<= $hlXn {
		qui replace x= hlWk[1,1] +   `k' * hlWk[1,3] in `L1'/`L2'

		projec in `L1'/`L2' 
        show	  `L1' `L2'

		loc L1=`L1'-1
		loc L2=`L2'-1
		if `L1'> 0  { qui replace y = y[_n+1]  in `L1'/`L2'}
		loc k =`k' +1
	}
end
*-------------
prog def show
*
	tempvar ylo yhi

	loc  L1 `1'
	loc  L2 `2'
	if !$hlPen {
		loc u=int(1+8*uniform())
		gph pen `u'
	}
	loc in "in `L1'/`L2'"
	qui { 
		g `ylo'= (gy <= lo)        `in'
		g `yhi'= (gy >= hi)        `in'
		replace isVi=`ylo'+2*`yhi' `in'
		capt assert isVi > 0       `in'
	}
	if _rc==0 { gph  vline gy gx in `L1'/`L2' }
	else 	  { 
		loc k `L1'
		loc L `L1'
		loc visi= isVi[`k']>0
		while `k'<`L2' {
			if `visi' {
				while `k'<`L2' &  isVi[`k'+1] { 
					loc k=`k'+1 
				}
				gph  vline gy gx in `L'/`k'
				
				if    `k'<`L2' {
					loc r=`k'+1
					cross `k' `r' 
					loc visi 0
				}
			}
				while `k'<`L2' & !isVi[`k'+1] { 
					loc k=`k'+1
					loc L=`k' 
				}
				if    `k'<`L2' {
					loc r=`k'+1
					cross `r' `k' 
					loc visi 1
					loc L=`r'			
				}
		}
	}			
	qui {
		replace lo= gy if `ylo' `in'
		replace hi= gy if `yhi' `in'
	}
end
*-------------
prog def cross
*
	loc  i  `1'
	loc  j  `2'

	loc low= isVi[`i']==1
	if !$hlPen { gph pen 7 }

	loc w0 = cond(`low', lo[`i'],hi[`i'])
	loc y0 = gy[`i']
	loc x0 = gx[`i']

	loc dw = cond(`low', lo[`j'],hi[`j']) - `w0'
	loc dy = gy[`j'] - `y0'
	loc dx = gx[`j'] - `x0'

	loc lam= - (`y0'-`w0')/(`dy'-`dw')
	loc cw = round(`w0' + `lam'*`dw', 1)
 
	loc cy = round(`y0' + `lam'*`dy', 1)
	loc cx = round(`x0' + `lam'*`dx', 1)

	gph line `y0' `x0' `cy' `cx'
end
*---------------------------------------------------
exit