pro get_amrv,t,err_t,lg,err_lg,feh,err_feh,logt_in,logg_in,logr_in,vmag_in,a,m,rad,Mv,boundarytouch=boundarytouch,noplot=noplot

;+
;	Using a grid from Bertelli's isochrones, we input  observed 
;	Teff, gravity, and [Fe/H] [t= Teff, lg= alog10(g); cgs units]
;	and we find out which values of a (age,  m (mass),  rad (radius)
;	and vmag (absolute magnitude)
;	are compatible for a star with the observed values.
;
;	IN: t	  - float	Teff (K)
;	    err_t - float	Uncertainty in Teff (K)
;	    lg	  - float	logg (g in cm -s^2)
;	    err_lg- float	Uncertainty in logg
;	    feh	  - float	[Fe/H]	= alog10(NFe/NH)-alog10(NFe/NH)_sun
;	   err_feh- float	Uncertainty in [Fe/H] (must be at least 0.25 dex)
;	   logt_in- fltarr(17,73,1800)	log10(Teff) as a function of [Fe/H],
;					log(age[ayr]), Mass/Msun (see below)
;	   logg_in- fltarr(17,73,1800)	log10(g) as a function of [Fe/H],
;					log(age[ayr]), Mass/Msun (see below)
;	   logr_in- fltarr(17,73,1800)	log10(radius) as a function of [Fe/H],
;					log(age[ayr]), Mass/Msun (see below)
;	   vmag_in- fltarr(17,73,1800)	Mv as a function of [Fe/H],
;					log(age[ayr]), Mass/Msun (see below)
;
;	OUT: a	- fltarr(3)	three-elements vector with an estimate
;			of log(Age[yr]), flanked by a lower an upper
;			limits ([lowerage,age,upperage]).
;	     m	- fltarr(3)	three-elements vector with an estimate
;			of the mass (initial - at birth), flanked by a lower 
;			an upper limits ([lowermass,mass,uppermass]).
;	     rad- fltarr(3)	three-elements vector with an estimate
;			of log10(radius), flanked by a lower an upper
;			limits ([lowerrad,rad,upperrad]).
;
;	     Mv-flatrr(3) 	three-elements vector with an estimate
;			of Mv , flanked by a lower an upper
;			limits ([lowerMv,Mv,upperMv]).
;
;
;	KEYWORDS: boundarytouch	This keyword returns a non-zero value when 
;		retrieved parameters gets solutions that touch the grid 
;		boundaries - an indication that the solution is likely biased.
;	
;
;	NOTE: logt and loggg are the grids read from the file 'bert.xdr'
;	computed with the routine mkgrid.pro (restore,'bert.xdr')
;
;	Example, the Sun: IDL> restore,'bert.xdr'
;			  IDL> get_amrv,5770.,150,4.44,0.1,0.0,0.25,logt,logg,logr,vmag,a,m,r
;			  IDL>print,a
;      -100000.      7.60000      9.29851	-> so, all we can say it's 
;						 likely younger than 2.0 Gyr
;			  IDL>print,m
;     0.701428      1.08546      1.46949	-> so it's M=1.1 +/- 0.4 Msun
;
; 	Carlos Allende Prieto, Sep/Dec 2002
;-
npar = n_params()
if (npar eq 0) then begin
	print,'get_amrv,t,err_t,lg,err_lg,feh,err_feh,LOGT,LOGG,LOGR,VMAG,a,m,r,Mv'
	return
endif 

logt=logt_in & logg=logg_in & logr=logr_in & vmag=vmag_in

lterr=float(err_t)/float(t)*alog10(exp(1)) 	; errorbars
lgerr=float(err_lg)
feherr=float(err_feh)
if (err_feh lt 0.25) then begin
	feherr=0.25
	print,'% GET_RAM: err_feh set to 0.25 dex in order to get mininum sampling'
endif


m_uniform=findgen(600)*0.006+0.6
la_uniform=findgen(73)*0.05+6.6
lz_uniform=findgen(17)*0.125-1.65

metallica=where(abs(lz_uniform-feh) le feherr)
if (max(metallica) eq -1) then begin
	metallica=where(abs(lz_uniform-feh) eq min(abs(lz_uniform-feh)))
endif
logt=logt(metallica,*,*)
logg=logg(metallica,*,*)
logr=logr(metallica,*,*)
vmag=vmag(metallica,*,*)

w=where((logt-alog10(t))^2/(3*lterr)^2+(logg-lg)^2/(3*lgerr)^2 le 1.0)
if (n_elements(w) lt 10) then begin
	a=[-1000,-1000,-1000]
	m=[-1000,-1000,-1000]
	rad=[-1000,-1000,-1000]
	mv=[-1000,-1000,-1000]
	print,'% GET_AM: Insufficient (< 10) solutions found'
	return
endif else begin

if not keyword_set(noplot) then begin
window,0
!p.multi=[0,0,0]
plot,logt(0,0,*),logg(0,0,*),xr=[4.4,3.5],yr=[5.0,1.0],xstyl=1,ystyl=1,/nodata,$
	xtitle='log Teff (K)',ytitle='logg (cm s-2)',charsi=1.5
	
print,'% GET_ZAM:',n_elements(w),' solutions found'
plotsym,3,3,/fill & loadct,12
plots,alog10(t),lg,psy=8,col=40
endif

aco=0 & zco=0 & mco=0 & peso=0
for i=0l,n_elements(w)-1 do begin
	donde,logt,w(i),coor
	zco=[zco,coor(0)] & aco=[aco,coor(1)] & mco=[mco,coor(2)]
	peso=[peso,$
	;(lterr/(logt(w(i))-alog10(t)))^2+(lgerr/(logg(w(i))-lg))^2]
	exp(-((logt(w(i))-alog10(t))/sqrt(2.)/lterr)^2 - $
		((logg(w(i))-lg)/sqrt(2.)/lgerr)^2)]
	if not keyword_set(noplot) then begin 
	plots,logt(coor(0),coor(1),coor(2)),logg(coor(0),coor(1),coor(2)),psy=3
	endif
	;print,n_elements(w),i
endfor

z=lz_uniform(zco(1:n_elements(zco)-1))
a=la_uniform(aco(1:n_elements(aco)-1))
m=m_uniform(mco(1:n_elements(mco)-1))
rad=logr(w)
Mv=vmag(w)
peso=peso(1:n_elements(peso)-1)

; normalizamos los pesos
peso=peso/total(peso)

if not keyword_set(noplot) then begin
	window,6
	!p.multi=[0,2,2] 
	ppl=1
endif else ppl=0

estimator_simple,a,peso,e,lower,upper,plot=ppl & a=[lower,e,upper]
;stop
estimator_simple,m,peso,e,lower,upper,plot=ppl & m=[lower,e,upper]
;stop
estimator_simple,rad,peso,e,lower,upper,plot=ppl & rad=[lower,e,upper]
;stop
estimator_simple,Mv,peso,e,lower,upper,plot=ppl & Mv=[lower,e,upper]
;stop

boundarytouch=0

if (max(where(z eq max(lz_uniform) or z eq min(lz_uniform))) gt -1) then begin
	print,'% GET_AMR: the retrieved Zs touch the grid boundaries'
	boundarytouch=1
endif

if (max(where(a eq max(la_uniform) or a eq min(la_uniform))) gt -1) then begin
	print,'% GET_AMR: the retrieved Ages touch the grid boundaries'
	boundarytouch=2	
endif

if (max(where(m eq max(m_uniform) or m eq min(m_uniform))) gt -1) then begin
	print,'% GET_AMR: the retrieved Masses touch the grid boundaries'
	boundarytouch=3	
endif

endelse

end