*! version 1.2.0 , apr 98, Guy van Melle   STB-45 gr30
*!
*! altitude plot: nearly a contour plot
*!
*! syntax: altitude  func  
*!		[,  x(lo hi step)  y(lo hi step) 
*!		Neg  Symb(int 4)  NQuant(int 16)  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 altitude
*----------------
	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 a 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
	}
	loc xmx= 32000 - 3
  	loc ymx= 23060 - 1500 * ($hlTex)
	init `xmx' `ymx' 

	gph open `sav'  
	 alti 
	 if $hlTex { hltex }
	gph close 
end
*----------------
prog def Get_info
*
	#delimit ;
	loc options "Xinfo(string) Yinfo(string) NQuant(int 16) Symb(int 4)
		Neg XMar(int 10) YMar(int 10) Text TMag(int 100) Penc
		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 "`xinfo'"=="" | "`yinfo'"=="" {
		di in red "xinfo, yinfo must  both  be provided"
		di in blu "x(lo,hi,step)  y(lo,hi,step) "
		exit 198
	}
	glo hlQua=   min( `nquant',20)
	glo hlSym=	  `symb'
	glo hlNeg=	("`neg'" !="")
	glo hlTex=	("`text'"!="")
	glo hlMag=	  `tmag'
	glo hlYma=	  `ymar'
	glo hlXma=	  `xmar'
	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 init
*
	loc xmx	= `1'
	loc ymx	= `2'
	mat hlWk= ($hlX \ $hlY )

	loc a = hlWk[1,2]-hlWk[1,1]
	loc da= hlWk[1,3]
	loc b = hlWk[2,2]-hlWk[2,1]
	loc db= hlWk[2,3]

	loc ok 0
	while !`ok' {
		glo hlXn= int(abs(`a'/`da'+0.1))
		glo hlYn= int(abs(`b'/`db'+0.1))
		glo hln = $hlXn + $hlYn + 1
		loc  ok = $hln <= $hlNmax
		if !`ok' {
			loc da=`da'* $hln / $hlNmax
			loc db=`db'* $hln / $hlNmax
		}
	}
	mat hlWk[1,3]=`da'
	mat hlWk[1,2]= hlWk[1,1]+$hlXn*`da'
	mat hlWk[2,3]=`db'
	mat hlWk[2,2]= hlWk[2,1]+$hlYn*`db'

* 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'
	}
	$hlprog x y z

*--- quietly set margins using bbox, apply scale and shift 

	qui {
		if $hlNeg { replace z= -z }
		loc yma= $hlYma/100*`ymx'
		loc xma= $hlXma/100*`xmx'
		loc ymb= `ymx'-`yma'
		loc xmb= `xmx'-`xma'

		set graphics off
		 gr y x, s(.) bbox(`yma', `xma', `ymb', `xmb', 800, 400, 0)
		set graphics on

		glo hlYfa= _result(5)
		glo hlYsh= _result(6)
		glo hlXfa= _result(7)
		glo hlXsh= _result(8)

		cap drop gx
		gen int  gx= int(x * $hlXfa + $hlXsh) 
		cap drop gy
		gen int  gy= int(y * $hlYfa + $hlYsh) 
	}
end
*------------
prog def alti
*
	tempvar  za zb ok
	tempname a

	loc  nq = $hlQua
	loc  sym= $hlSym
	qui {
		g  `za'=z
		g  `zb'=0
		g  `ok'=1
		loc nc = `nq'-1
		loc dc = 800/`nq'
		mat `a'= J(1,`nc',0)

		_pctile `za', n(`nq')
		loc k 0
		while `k'<`nc' {
			loc k=`k'+1
			mat `a'[1,`k']=_result(`k')
		}
			gr 
			replace `ok'= (z<= `a'[1,1])
			replace `za'= 20
			gphit  1  `za' `zb' `ok' 
		loc k 1
		replace `zb'= `sym'
		while `k'<`nc' {
			replace `ok'= (z > `a'[1,`k'])
			loc k=`k'+1
			replace `ok'= (z<= `a'[1,`k']) if `ok'
			replace `za'= round(`za'+`dc',1)
			gphit `k' `za' `zb' `ok'  
		}
			replace `ok'= (z > `a'[1,`k'])
			loc k=`k'+1
			replace `za'= round(`za'+`dc',1)
			gphit `k' `za' `zb' `ok'  
	}
end
*-------------
prog def gphit
*
	loc k = `1'
	loc za  `2'
	loc zb  `3'
	loc ok  `4'
	loc y1 = 1000+int((`k'-1)*20000/$hlQua)
	loc y2 = `y1'-200
	loc x1 = 30200
	loc x2 = `x1'+1200
	loc p = mod((`k'-1),8)+1
	loc z =`za'[1]

	gph  pen `p'	
	gph vpoint  gy   gx  `za' `zb' if `ok'
	gph  text  `y1' `x1' 0 0 `k'=
	gph  point `y2' `x2' `z'  $hlSym
end
*------------