********************************************************************************************************************************
* pbpoint fits the Preece and Baines 1978 family of growth curves.
* version 1.0
* Author: Adrian Sayers
* Date: 06.03.2013
*
********************************************************************************************************************************
	prog define pbpoint , sclass
		version 9.2
		syntax varlist(min=2 max=2 numeric) [if] [in] , [Model( integer 1) iterate(integer 100)  ]

tokenize `varlist'
	tempname outcome time
		gen `outcome' = `1'
		gen `time' = `2'

		if `model'==1 |`model'==2 {


		if `model'==1 {
		qui nlcom 	(h1: _b[/h1]) (htheta: _b[/htheta]) (s0: _b[/s0]) (s1: _b[/s1]) (theta: _b[/theta]) (g: 1) , post
						}
		if `model'==2 {
		qui nlcom 	(h1: _b[/h1]) (htheta: _b[/htheta]) (s0: _b[/s0]) (s1: _b[/s1]) (theta: _b[/theta]) (g: _b[/ga]) , post
					}


# delimit ;
noisily di _n "Calculating Point Estimates" _n ;

local nliterate = `iterate' ;
* Calculate values of s ;

capture nlcom 	(s_neg: ( - _b[g] * (_b[s0] + _b[s1]) - sqrt(_b[g]^2 * (_b[s0] - _b[s1])^2 - 4*_b[g] * _b[s0] * _b[s1] )) / ( - 2 - 2 * _b[g]))
		(s_pos: ( - _b[g] * (_b[s0] + _b[s1]) + sqrt(_b[g]^2 * (_b[s0] - _b[s1])^2 - 4 * _b[g] * _b[s0] * _b[s1] )) / ( - 2 - 2 * _b[g]))
		(h1: _b[h1]) (htheta: _b[htheta]) (s0: _b[s0]) (s1: _b[s1]) (theta: _b[theta]) (g: _b[g])		, post iterate(`nliterate');

if _rc!=0 {; di as err "Point Estimates Fail" ; exit _rc ; } ;

* Calculate age at phv , and take off ;

capture nlcom	(t_neg: _b[theta] + (log( - (_b[s_neg] - _b[s1]) / (_b[s_neg] - _b[s0])) / (_b[g] * (_b[s0] - _b[s1] ))))
(t_pos: _b[theta] + (log( - (_b[s_pos] - _b[s1]) / (_b[s_pos] - _b[s0])) / (_b[g] * (_b[s0] - _b[s1]))))
(h1: _b[h1]) (htheta: _b[htheta]) (s0: _b[s0]) (s1: _b[s1]) (theta: _b[theta]) (g: _b[g]) (s_neg: _b[s_neg]) (s_pos: _b[s_pos]) , post iterate(`nliterate') ;

if _rc!=0 {; di as err "Point Estimates Fail" ; exit _rc ; } ;

* Calculate Max velocity ;

capture nlcom 	(ht_vel_pos: _b[h1] - ((_b[h1] - _b[htheta]) / (((0.5 * exp(_b[g] * _b[s0] * (_b[t_pos]-_b[theta]))) + (0.5 * exp(_b[g] * _b[s1] * (_b[t_pos] - _b[theta]))))^(1 / _b[g]))))
		(ht_vel_neg: _b[h1] - ((_b[h1] - _b[htheta]) / (((0.5 * exp(_b[g] * _b[s0] * (_b[t_neg]-_b[theta]))) + (0.5 * exp(_b[g] * _b[s1] * (_b[t_neg] - _b[theta]))))^(1 / _b[g]))))
	(h1: _b[h1]) (htheta: _b[htheta]) (s0: _b[s0]) (s1: _b[s1]) (theta: _b[theta]) (g: _b[g])
	(s_neg: _b[s_neg]) (s_pos: _b[s_pos])  (t_neg: _b[t_neg]) (t_pos: _b[t_pos])
	, post iterate(`nliterate') ;

if _rc!=0 {; di as err "Point Estimates Fail" ; exit _rc ; } ;


capture nlcom 	(vel_pos:  _b[s_pos] * (_b[h1] - _b[ht_vel_pos]) )
		(vel_neg:  _b[s_neg] * (_b[h1] - _b[ht_vel_neg]) )
		(h1: _b[h1]) (htheta:_b[htheta]) (s0:_b[s0]) (s1:_b[s1]) (theta:_b[theta]) (g:_b[g])
		(s_neg: _b[s_neg]) (s_pos: _b[s_pos]) (t_neg: _b[t_neg]) (t_pos: _b[t_pos]) (ht_vel_pos: _b[ht_vel_pos]) (ht_vel_neg: _b[ht_vel_neg])
		, post iterate(`nliterate') ;

if _rc!=0 {; di as err "Point Estimates Fail" ; exit _rc ; } ;

local h1  = _b[h1] ;
local htheta	= _b[htheta] ;

if `h1' > `htheta' { ;
 capture noisily nlcom 	 (age_at_take_off: _b[t_pos]) 	(ht_at_take_off: _b[ht_vel_pos]) 	(vel_at_take_off: _b[vel_pos])
   				 (age_at_phv : _b[t_neg]) 		(ht_at_phv: _b[ht_vel_neg]) 		(phv: _b[vel_neg])
 				 (age_at_finish_off: 0 ) 		(ht_at_finish_off: 0 ) 				(vel_at_finish_off: 0 ), post  level(`level') iterate(`nliterate') ;
 					 if _rc!=0 {; di as err "Point Estimates Fail" ; exit _rc ; } ;
 				 } ;

if `h1' < `htheta' { ;
 capture noisily nlcom 	 (age_at_take_off: 0) 				(ht_at_take_off: 0) 				(vel_at_take_off: 0)
  				 (age_at_phv : _b[t_pos]) 			(ht_at_phv: _b[ht_vel_pos]) 		(phv: _b[vel_pos])
 				 (age_at_finish_off: _b[t_neg] ) 	(ht_at_finish_off: _b[ht_vel_neg] ) (vel_at_finish_off: _b[vel_neg] ), post  level(`level') iterate(`nliterate') ;
						if _rc!=0 {; di as err "Point Estimates Fail" ; exit _rc ; } ;
 				  				 		} ;
estimates store pb_point ;

#delimit cr
}

end