*! version 1.2 am 8Jul98              [STB-51: gr40]
program define contour
    version 5.0
    local varlist "req ex min(2) max(100)"
    local options "SAVING(string) LTitle(string) BTitle(string) TICKS(integer 5) TItle(string) Box(string) CONtour(string) PEN(string) TEXT LEGEND"
    parse "`*'"
    parse "`varlist'", parse(" ")

* the default colour is white gphpen 3
global def_col=3

di "$S_TIME"

macno "`varlist'" " "
global XDIM=$macn
global YDIM=_N

di "$XDIM columns by $YDIM"

*check variables

parse "`varlist'", parse(" ")
local i 1
while `i'<=$macn {
	if "``i''"=="x`i'" {
	}
	else {
		di "You must have the x variables incrementing by one ``i'' should be x`i'"
		di "OR you must call each variable x(integer) not ``i''"
		exit(666)		
	}

local i=`i'+1
}



/* Set up graph box variables

	leftx			rightx
topy	-------------------------
	|			|
	|			|
height	|			|
	|			|
	|			|
	|			|
boty	-------------------------
	<-	width		->
*/

global max_by=23063
global max_rx=32000

global topy=2000
global boty=21000
global leftx=2500
global rightx=30000
global width=$rightx-$leftx
global height=$boty-$topy

if "`box'" ~= "" {
	parse "`box'", parse(",")
	if(`1'<2000) {
		di "WARNING Bounding box wrong top y `1'<2000"
		exit(666)
	}
	if(`3'>21000) {
		di "WARNING Bounding box wrong bottom y `3'>21000"
		exit(666)
	}
	if(`5'<2500) {
		di "WARNING Bounding box wrong left x `5'<2500"
		exit(666)
	}
	if(`7'>30000) {
		di "WARNING Bounding box wrong right x `7'>30000"
		exit(666)
	}
	if(`1'>=`3' | `5'>=`7') {
		di "WARNING Bounding box wrong Is `1'>=`3'? OR is `5'>=`7'?"
		exit(666)
	}
	else {
		global topy=`1'
		global boty=`3'
		global leftx=`5'
		global rightx=`7'
		global width=$rightx-$leftx
		global height=$boty-$topy
	}

}

global dleftx=2000
global drightx=30000
global dtopy=2000
global dboty=20000

global xwid= int(($drightx-$dleftx)/$XDIM)
global ywid= int(($dboty-$dtopy)/$YDIM)

*
* Saving the file or not.
* NB if file exists it is deleted!!
*

	if ("`saving'"~="") {
		cap confirm new file `saving'.gph
		if _rc~=0 {
			di "Deleteing file...`saving'.gph"
			!rm `saving'.gph
		}
		di "Saving file... `saving'.gph"
		gph open, saving(`saving')
	}
	else { gph open }


if  "`contour'"~="" {
	local mpen = 1
	parse "`contour'", parse(",")
	while "``mpen''" ~="" {
		parse "`pen'", parse(",")
		if "``mpen''"~="" {
			gph pen ``mpen''
		}
		else { gph pen 2 }

		parse "`contour'", parse(",")
		draw_con ``mpen''
		if "`legend'" ~= "" {
			local xxx = $drightx+100
			local yyy = $boty-`mpen'*350+2000
			gph text `yyy' `xxx'  0 0 ``mpen''
		}
		local mpen = `mpen'+2
	}
}
else {
preserve
stack x*, into(ade) clear
qui summ ade,de
local con1=_result(9)
local con2=_result(10)
local con3=_result(11)
restore

local pen 2
local ii=1
while `ii'<4 {

gph pen `pen'
draw_con `con`ii''
local pen=`pen'+2

	if "`legend'" ~= "" {
		local xxx = $drightx+1000
		local yyy = $boty-`pen'*350+3000
		gph text `yyy' `xxx'  0 0 `con`ii''
	}
local ii=`ii'+1
}
}

gph pen $def_col

draw_axe $dleftx $drightx $dtopy $dboty

draw_tic `ticks'

	local xxx = ($drightx-$dleftx)/2+$dleftx
	local yyy = $dboty+1750
	gph text `yyy' `xxx' 0 0 `btitle'
	local yyy = $dtopy -750
	gph text `yyy' `xxx' 0 0 `title'
	local xxx= $dleftx
	gph text `yyy' `xxx' 0 0 `ltitle'

if "`text'"~="" {
	tempvar yy1 xx1 txt
	gen `yy1'=0
	gen `xx1'=0
	gen str10 `txt'=""
 
	local i=1
	while `i'<=$XDIM {

		qui replace `yy1'=cond(1,int($dtopy+_n*$ywid - $ywid/2)  ,0)
		qui replace `xx1'=cond(1,int($dleftx+`i'*$xwid - $xwid/2) ,0)
		qui replace `txt'= string(x`i'[_n])
		gph vtext `yy1' `xx1' `txt'
	local i=`i'+1
	}
}

gph close
	

if "`if'"~="" {
	restore
}
di "$S_TIME"
end

program define draw_axe

gph line `3' `1' `3' `2'
gph line `4' `1' `4' `2'
gph line `3' `2' `4' `2'
gph line `3' `1' `4' `1'

end

program define draw_tic
tempvar xtic1 ytic1 xtic2 ytic2  txt

gen `xtic1'=0
gen `xtic2'=0
gen `ytic1'=0
gen `ytic2'=0
gen str4 `txt'=""

local i=1
while `i'<=$XDIM {

	if mod(`i',`1')==0 {
	local xtic1 = int($dleftx+(`i'-1)*$xwid) 
	local xtic2 = int($dleftx+(`i')*$xwid) 
	local ytic1 = int($dboty)  
	local ytic2 = int($dboty+300) 
	gph line `ytic1' `xtic1' `ytic2' `xtic1' 
	gph line `ytic1' `xtic2' `ytic2' `xtic2'
	local ytic1 = int($dboty) +900	 
	local xtic2 = int($dleftx+(`i'-0.5)*$xwid) 
	gph text `ytic1' `xtic2' 0 0 `i'
	}
local i = `i'+1
}

local i=1
while `i'<=$YDIM {

	if mod(`i',`1')==0 {
	local ytic1 = int($dtopy+(`i'-1)*$ywid) 
	local ytic2 = int($dtopy+(`i')*$ywid) 
	local xtic1 = int($dleftx)  
	local xtic2 = int($dleftx-300) 
	gph line `ytic1' `xtic1' `ytic1' `xtic2' 
	gph line `ytic2' `xtic1' `ytic2' `xtic2'
	local xtic2 = int($dleftx-300) -700
	local ytic1 = int($dtopy+(`i'-0.25)*$ywid) 
	gph text `ytic1' `xtic2' 0 0 `i'
	}
local i = `i'+1
}

end

*
* Drawing a contour taking one value and drawing all the edges
* that are needed to draw continuous lines.
*
*

program define draw_con

tempvar xline yline yy1 yy2 yy3 xx1 xx2 xx3 txt

gen `yy1'=0
gen `yy2'=0
gen `yy3'=0
gen `xx1'=0
gen `xx2'=0
gen `xx3'=0
gen str4 `txt'=""
local i 1
while `i' <= $XDIM {
	tempvar x`i'
	qui gen `x`i''=x`i'>`1'
local i=`i'+1
}

mat temp=J(10,10,0)

local i=1
while `i'<=$XDIM {

/*TEXT
		qui replace `yy1'=cond(1,int($dtopy+_n*$ywid - $ywid/2)  ,0)
		qui replace `xx1'=cond(1,int($dleftx+`i'*$xwid - $xwid/2) ,0)
		qui replace `txt'= string(`x`i''[_n]) 
		gph vtext `yy1' `xx1' `txt'
*/
		local il=`i'-1
		local ir = `i'+1

*Vertical lines1
		qui replace `yy1'=cond(`x`i''[_n]>`x`ir''[_n] & `i'<$XDIM , /*
			*/ cond(`x`ir''[_n-1]==1,  int($dtopy+_n*$ywid - $ywid/2) ,/*
			*/ cond(`x`i''[_n-1]==1 | _n==1, int($dtopy+_n*$ywid-$ywid/2)- $ywid/2, int($dtopy+_n*$ywid-$ywid/2)) /*
			*/ ) /*
			*/ ,0)
		qui replace `yy2'=cond(`x`i''[_n]>`x`ir''[_n] ,/*
			*/ cond(`x`ir''[_n+1]==1 ,  int($dtopy+_n*$ywid - $ywid/2) ,/*
			*/ cond(`x`i''[_n+1]==1 | _n==_N, int($dtopy+_n*$ywid-$ywid/2)+$ywid/2, int($dtopy+_n*$ywid-$ywid/2)) /*
			*/ ) /*
			*/ , 0)

		qui replace `xx1'=cond(`x`i''[_n]>`x`ir''[_n] ,int($dleftx+`i'*$xwid - $xwid/2)+$xwid/2,0)
		qui replace `xx2'=cond(`x`i''[_n]>`x`ir''[_n] ,int($dleftx+`i'*$xwid - $xwid/2)+$xwid/2 ,0)
		gph vpoly `yy1' `xx1' `yy2' `xx2' if `yy1'~=0

*Vertical lines2
		qui replace `yy1'=cond(`x`i''[_n]<`x`ir''[_n] & `i'<$XDIM , /*
			*/ cond(`x`ir''[_n-1]==0 ,  int($dtopy+_n*$ywid - $ywid/2) ,/*
			*/ cond(`x`i''[_n-1]==0  | _n==1, int($dtopy+_n*$ywid-$ywid/2)- $ywid/2, int($dtopy+_n*$ywid-$ywid/2)) /*
			*/ ) /*
			*/ ,0)
		qui replace `yy2'=cond(`x`i''[_n]<`x`ir''[_n] ,/*
			*/ cond(`x`ir''[_n+1]==0 ,  int($dtopy+_n*$ywid - $ywid/2) ,/*
			*/ cond(`x`i''[_n+1]==0 | _n==_N, int($dtopy+_n*$ywid-$ywid/2)+$ywid/2, int($dtopy+_n*$ywid-$ywid/2) ) /*
			*/ ) /*
			*/ , 0)

		qui replace `xx1'=cond(`x`i''[_n]<`x`ir''[_n] ,int($dleftx+`i'*$xwid - $xwid/2)+$xwid/2,0)
		qui replace `xx2'=cond(`x`i''[_n]<`x`ir''[_n] ,int($dleftx+`i'*$xwid - $xwid/2)+$xwid/2 ,0)
		gph vpoly `yy1' `xx1' `yy2' `xx2' if `yy1'~=0

*horiz1

		qui replace `yy1'=cond(`x`i''[_n]<`x`i''[_n+1] &_n<_N, /*
			*/ cond(`x`il''[_n]==1,  int($dtopy+_n*$ywid - $ywid/2) ,/*
			*/ cond(`x`il''[_n+1]==1, int($dtopy+_n*$ywid-$ywid/2)+$ywid/2, int($dtopy+_n*$ywid-$ywid/2)+$ywid) /*
			*/ ) /*
			*/ ,0)
		qui replace `yy2'=cond(`x`i''[_n]<`x`i''[_n+1] &_n<_N,int($dtopy+_n*$ywid - $ywid/2) +$ywid/2,0)
		qui replace `yy3'=cond(`x`i''[_n]<`x`i''[_n+1] &_n<_N,/*
			*/ cond(`x`ir''[_n]==1,  int($dtopy+_n*$ywid - $ywid/2) ,/*
			*/ cond(`x`ir''[_n+1]==1, int($dtopy+_n*$ywid-$ywid/2)+$ywid/2, int($dtopy+_n*$ywid-$ywid/2)+$ywid) /*
			*/ ) /*
			*/ ,0)

		qui replace `xx1'=cond(`x`i''[_n]<`x`i''[_n+1] &_n<_N,int($dleftx+`i'*$xwid - $xwid/2)-$xwid/2,0)
		qui replace `xx2'=cond(`x`i''[_n]<`x`i''[_n+1] &_n<_N,int($dleftx+`i'*$xwid - $xwid/2)  ,0)
		qui replace `xx3'=cond(`x`i''[_n]<`x`i''[_n+1] &_n<_N,int($dleftx+`i'*$xwid - $xwid/2) +$xwid/2,0)
		gph vpoly `yy1' `xx1' `yy2' `xx2' `yy3' `xx3' if `yy1'~=0


*horiz2

		qui replace `yy1'=cond(`x`i''[_n]>`x`i''[_n+1] &_n<_N, /*
			*/ cond(`x`il''[_n+1]==1,  int($dtopy+_n*$ywid - $ywid/2)+$ywid ,/*
			*/ cond(`x`il''[_n]==1, int($dtopy+_n*$ywid-$ywid/2)+$ywid/2, int($dtopy+_n*$ywid-$ywid/2)) /*
			*/ ) /*
			*/ ,0)
		qui replace `yy2'=cond(`x`i''[_n]>`x`i''[_n+1] &_n<_N,int($dtopy+_n*$ywid - $ywid/2) +$ywid/2,0)
		qui replace `yy3'=cond(`x`i''[_n]>`x`i''[_n+1] &_n<_N,/*
			*/ cond(`x`ir''[_n+1]==1,  int($dtopy+_n*$ywid - $ywid/2)+$ywid ,/*
			*/ cond(`x`ir''[_n]==1, int($dtopy+_n*$ywid-$ywid/2)+$ywid/2, int($dtopy+_n*$ywid-$ywid/2)) /*
			*/ ) /*
			*/ ,0)

		qui replace `xx1'=cond(`x`i''[_n]>`x`i''[_n+1] &_n<_N,int($dleftx+`i'*$xwid - $xwid/2)-$xwid/2,0)
		qui replace `xx2'=cond(`x`i''[_n]>`x`i''[_n+1] &_n<_N,int($dleftx+`i'*$xwid - $xwid/2)  ,0)
		qui replace `xx3'=cond(`x`i''[_n]>`x`i''[_n+1] &_n<_N,int($dleftx+`i'*$xwid - $xwid/2) +$xwid/2,0)
		gph vpoly `yy1' `xx1' `yy2' `xx2' `yy3' `xx3' if `yy1'~=0


local i=`i'+1
}

end

program define macno

global macn=0
global last=0
local tmp_str = ""
local tmp_str "`2'"
parse "`1'", parse("`2'")

while "`1'"~="" {
	while "`1'"=="`tmp_str'" {
		mac shift
		global last = $last+1
	}
	global macn = $macn+1
	mac shift
global last= $last+1

}

end