#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#define OBS_TYPE_MASS 1
#define OBS_TYPE_RADIUS 2
#define OBS_TYPE_DENSITY 3
#define OBS_TYPE_LOGDENSITY 4
#define OBS_TYPE_TEFF 5
#define OBS_TYPE_VMAG 6
#define OBS_TYPE_JMAG 7
#define OBS_TYPE_HMAG 8
#define OBS_TYPE_KMAG 9
#define OBS_TYPE_ABSVMAG 10
#define OBS_TYPE_ABSJMAG 11
#define OBS_TYPE_ABSHMAG 12
#define OBS_TYPE_ABSKMAG 13
#define OBS_TYPE_VJCOLOR 14
#define OBS_TYPE_VHCOLOR 15
#define OBS_TYPE_VKCOLOR 16
#define OBS_TYPE_JHCOLOR 17
#define OBS_TYPE_JKCOLOR 18
#define OBS_TYPE_HKCOLOR 19
#define OBS_TYPE_DISTMOD 20
#define OBS_TYPE_DIST 21
#define OBS_TYPE_AV 22
#define OBS_TYPE_EBV 23
#define OBS_TYPE_LOGG 24
#define OBS_TYPE_LOGGSOLAR 25
#define OBS_TYPE_DENSITYSOLAR 26
#define OBS_TYPE_LOGDENSITYSOLAR 27
#define OBS_TYPE_PARALLAX 28
#define OBS_TYPE_ANGDIAM 29

#define NGR_PASS_STOP_LIMIT 2

/* Parameters from reddenning law */
#define REDLAW_AB 1.36
#define REDLAW_AJ 0.282
#define REDLAW_AH 0.179
#define REDLAW_AK 0.114

#define LOGGSUN 4.438  /* CGS units */

#define SOLARDENSITY 1.4115  /* CGS units */

#define MIN_LRHO_LIMIT 0  /* Minimum log10(rho) which the model will consider */
#define MAX_LRHO_LIMIT 2.5 /* Maximum log10(rho) which the model will consider */

#define MAX_INPUT_MASS 0.8
#define MAX_INPUT_RADIUS 1.0

#define BADCHI2 -99999.0

double ModelMedian[22] = { 0.054, -0.2989, 0.0815, -0.18280, 32.6, -7.09, 1.24, -7.690, 570, 16.5, 0.111, 100, -5.4, -0.1105, -130, -3.46, -0.1615, -80, -4.83, -0.1757, 0.0091, 0.169 };

double ModelStddev[22] = { 0.049, 0.0084, 0.0020, 0.00075, 3.3, 0.44, 0.11, 0.034, 380, 2.9, 0.013, 190, 1.2, 0.0046, 150, 1.13, 0.0044, 140, 0.99, 0.0039, 0.0012, 0.025 };

typedef struct {
  int Nobs;
  double *obsval;
  double *errval;
  int *obstype;
  char *isupper;
  char *islower;
  double *upperlim;
  double *lowerlim;
} _ObsData;

typedef struct {
  double logrho;
  double dlR;
  double dMbol;
  double lM;
  double lR;
  double lTeff;
  double Mbol;
  double BCV;
  double BCJ;
  double BCH;
  double BCK;
  double absVmag;
  double absJmag;
  double absHmag;
  double absKmag;
  double AV;
  double distmod;
} _ModelData;

double GetLnLikelihood(_ObsData *data, double *paramvec, _ModelData *m, double isAV, double isdist, double maxAV);

#define TINYERRTERM 1.0e-5

#define NEMPIRICALMODELPARAMS 22
#define INDEX_a0 0
#define INDEX_a1 1
#define INDEX_a2 2
#define INDEX_a3 3
#define INDEX_b0 4
#define INDEX_b1 5
#define INDEX_b2 6
#define INDEX_b3 7
#define INDEX_aBCV0 8
#define INDEX_aBCV1 9
#define INDEX_aBCV2 10
#define INDEX_aBCJ0 11
#define INDEX_aBCJ1 12
#define INDEX_aBCJ2 13
#define INDEX_aBCH0 14
#define INDEX_aBCH1 15
#define INDEX_aBCH2 16
#define INDEX_aBCK0 17
#define INDEX_aBCK1 18
#define INDEX_aBCK2 19
#define INDEX_Sr 20
#define INDEX_Sm 21
typedef struct {
  double v[NEMPIRICALMODELPARAMS];
} _EmpiricalModelParams;

#ifdef SWAP
#undef SWAP
#endif
#define SWAP(a,b) temp=(a);(a)=(b);(b)=temp

#define DEF_QUICKSELECT(type) \
type quickselect(int k, int n, type *arr) \
{ \
  int i, ir, j,m, mid; \
  type a, temp; \
  k--; \
  m = 0; \
  ir = n-1; \
  for(;;) { \
    if(ir<=m+1) { \
      if(ir == m+1 && arr[ir] < arr[m]) { \
	SWAP(arr[m],arr[ir]); \
      } \
      return arr[k]; \
    } else { \
      mid=(m+ir) >> 1; \
      SWAP(arr[mid],arr[m+1]); \
      if(arr[m] > arr[ir]) { \
	SWAP(arr[m],arr[ir]); \
      } \
      if(arr[m+1] > arr[ir]) { \
	SWAP(arr[m+1],arr[ir]); \
      } \
      if(arr[m] > arr[m+1]) { \
	SWAP(arr[m],arr[m+1]); \
      } \
      i=m+1; \
      j=ir; \
      a=arr[m+1]; \
      for(;;) { \
	do i++; while(arr[i] < a); \
	do j--; while(arr[j] > a); \
	if(j < i) break; \
	SWAP(arr[i],arr[j]); \
      } \
      arr[m+1]=arr[j]; \
      arr[j]=a; \
      if(j >= k) ir=j-1; \
      if(j <=k) m=i; \
    } \
  } \
}

DEF_QUICKSELECT(double);

#undef SWAP

/* Computes the median of a vector */
double ComputeMedian(double *data, int n)
{
  double temp;
  if(n == 1)
    return data[0];
  if(n % 2 == 0) {
    temp = quickselect((n/2),n,data);
    temp += quickselect((n/2)+1,n,data);
    temp = temp/2.0;
  }
  else {
    temp = quickselect((n/2)+1,n,data);
  }
  return temp;
}

/* Computes the standard deviation of a vector */
double ComputeStddev(double *data, int n) 
{
  int i;
  double ave = 0.;
  double ave2, term, stddevout;
  /* First get the average */
  for(i=0; i < n; i++) {
    ave += data[i]/n;
  }
  /* Now get the standard deviation */
  ave2 = 0.;
  for(i=0; i < n; i++) {
    term = (data[i] - ave);
    ave2 += term*term;
  }
  stddevout = sqrt(ave2/((double)(n - 1)));
  return stddevout;
}

#ifdef TINY
#undef TINY
#endif
#define TINY 1.0e-10
#define VARTOOLS_AMOEBA_DEFAULT_MAXSTEPS 5000

#ifdef SWAP
#undef SWAP
#endif

#define SWAP(a,b) swap=(a);(a)=(b);(b)=swap

int amoeba(double **p, double *y, int *ia, int ndim, double ftol, int *nfunk, int maxeval, _ObsData *data, _ModelData *m, double isAV, double isdist, double maxAV)
/* This function runs the amoeba downhill simplex algorithm.
   p - a 2-d array storing the simplex
   y - a vector storing the chi2 values for each row in the simplex
   ia - an integer array storing flags that indicate whether or
        not the parameter is to be varied. The function will allocate memory for
        the ia array, and update it, for each parameter added. This array will
        be passed to amoeba.
   ndim - the number of parameters in the simplex.
   ftol - the fractional in chi2 to decide when to stop the minimization.
   funk - a pointer to a function which provides chi2.
      funk has the expected format:
      double funk(double *param, int ndim, int N, double *t, double *mag,
                  double *sig, void *userparams)
           where param is a set of parameter values to calculate chi2 for,
           and the function should return chi2.
      The other input parameters to funk are as described below.
   nfunk - on return, this stores the number of function calls executed by
           amoeba.
   maxeval - maximum number of function evaluations to perform before giving
             up. If this is <= 0 then the default number is used.
   The following terms are passed directly to funk and not actually used
   elsewhere within amoeba, they can, but do not have to, take on the context
   given below:
    N - number of points in the time series that is being fit.
    t - times of observation
    mag - set of observed magnitude values being fit
    sig - set of magnitude uncertainties
    userparams - structure containing other user-defined parameters.
         must be cast to the appropriate type by funk.
*/
{
  double amotry(double **p, double *y, int *ia, double *psum, int ndim, int ndimused, int ihi, double fac, _ObsData *data, _ModelData *m, double isAV, double isdist, double maxAV);
  int i, ihi, ilo, inhi, j, mpts, ndimused;
  double rtol, sum, swap, ysave, ytry, *psum;

  if(maxeval <= 0) maxeval = VARTOOLS_AMOEBA_DEFAULT_MAXSTEPS;

  ndimused = 0;
  for(j=0;j<ndim;j++)
    if(ia[j])
      ndimused++;

  if(!ndimused)
    return 1;

  mpts = ndimused+1;

  if((psum = (double *) malloc(ndim * sizeof(double))) == NULL)
    {
      fprintf(stderr,"Memory allocation error\n");
      exit(1);
    }
  *nfunk = 0;
  for (j=0;j<ndim;j++) {
    for(sum=0.0,i=0;i<mpts;i++) sum += p[i][j];
    psum[j]=sum;
  }
  for(;;) {
    ilo = 1;
    ihi = y[0] > y[1] ? (inhi=1,0) : (inhi=0,1);
    for (i=0;i<mpts;i++) {
      if(y[i] <= y[ilo]) ilo = i;
      if(y[i] > y[ihi]) {
	inhi = ihi;
	ihi = i;
      } else if(y[i] > y[inhi] && i != ihi) inhi = i;
    }
    rtol = 2.0*fabs(y[ihi]-y[ilo])/(fabs(y[ihi])+fabs(y[ilo])+TINY);
    if(rtol < ftol) {
      SWAP(y[0],y[ilo]);
      for(i=0;i<ndim;i++) {SWAP(p[0][i],p[ilo][i]);}
      break;
    }
    if (*nfunk >= maxeval)
      return 1;
    //error(ERR_TOOMANYAMOEBAITERATIONS);

    *nfunk += 2;
    ytry = amotry(p,y,ia,psum,ndim,ndimused,ihi,-1.0, data, m, isAV, isdist, maxAV);
    
    if(!(ytry <= BADCHI2) && ((y[ilo] <= BADCHI2) || ytry <= y[ilo]))
      ytry = amotry(p,y,ia,psum,ndim,ndimused,ihi,2.0,data, m, isAV, isdist, maxAV);
    else if((ytry <= BADCHI2) || ytry >= y[inhi])
      {
	ysave=y[ihi];
	ytry = amotry(p,y,ia,psum,ndim,ndimused,ihi,0.5,data, m, isAV, isdist, maxAV);
	if ((ytry <= BADCHI2) || ytry >= ysave) {
	  for (i = 0; i<mpts;i++) {
	    if (i != ilo) {
	      for(j=0;j<ndim;j++)
		{
		  if(ia[j])
		    {
		      p[i][j] = psum[j] = 0.5*(p[i][j] + p[ilo][j]);
		    }
		  else
		    {
		      psum[j] = p[i][j];
		    }
		}
	      y[i]=GetLnLikelihood(data, psum, m, isAV, isdist, maxAV);
	    }
	  }
	  *nfunk += ndim;
	  for (j=0;j<ndim;j++) {
	    for(sum=0.0,i=0;i<mpts;i++) sum += p[i][j];
	    psum[j]=sum;
	  }
	}
      } else --(*nfunk);
  }
  free(psum);
  return 0;
}

double amotry(double **p, double *y, int *ia, double *psum, int ndim, int ndimused, int ihi, double fac, _ObsData *data, _ModelData *m, double isAV, double isdist, double maxAV)
{
  int j;
  double fac1, fac2, ytry, *ptry;
  if((ptry = (double *) malloc(ndim * sizeof(double))) == NULL)
    {
      fprintf(stderr,"Memory Allocation Error\n");
      exit(1);
    }
  fac1 = (1.0 - fac)/ndimused;
  fac2 = fac1 - fac;
  for(j=0;j<ndim;j++)
    {
      if(ia[j])
	ptry[j]=psum[j]*fac1 - p[ihi][j]*fac2;
      else
	ptry[j]=p[ihi][j];
    }
  ytry=GetLnLikelihood(data, ptry, m, isAV, isdist, maxAV);
  if (!(ytry <= BADCHI2) && (ytry < y[ihi])) {
    y[ihi]=ytry;
    for (j = 0; j< ndim;j++) {
      if(ia[j])
	{
	  psum[j] += ptry[j] - p[ihi][j];
	  p[ihi][j] = ptry[j];
	}
    }
  }
  free(ptry);
  return ytry;
}

void amoeba_initializesimplexchi2(int Nparameters, int Ntovary, double **p, double **chi2vals, _ObsData *data, _ModelData *m, double isAV, double isdist, double maxAV)
/* This function sets up the initial chi2 vector for the amoeba simplex.
   This is called after the simplex has been setup using 
   incrementparameters_foramoeba.
   Nparameters - the number of parameters in the simplex.
   Ntovary - the total number of parameters that are free to vary.
   p - A 2-d array storing the simplex.
   chi2vals - A pointer to a 1-d vector which will store the chi2 values for
     each point on the simplex. The program will allocate memory for this 
     vector.
   funk - a pointer to a function which provides chi2.
      funk has the expected format:
      double funk(double *param, int Nparameters, int N, double *t, double *mag,
                  double *sig, void *userparams)
           where param is a set of parameter values to calculate chi2 for,
           and the function should return chi2.
      The other input parameters to funk are as described below.
   The following terms are passed directly to funk and not actually used
   elsewhere within this function, they can, but do not have to, 
   take on the context given below:
    N - number of points in the time series that is being fit.
    t - times of observation
    mag - set of observed magnitude values being fit
    sig - set of magnitude uncertainties
    userparams - structure containing other user-defined parameters.
         must be cast to the appropriate type by funk.
*/
{
  int i, j;
  double gaussian(void);
  if(((*chi2vals) = (double *) malloc((Ntovary + 1)*sizeof(double))) == NULL)
    {
      fprintf(stderr,"Memory Allocation Error\n");
      exit(1);
    }
  for(i=0; i < Ntovary + 1; i++) {
    (*chi2vals)[i] = GetLnLikelihood(data, p[i], m, isAV, isdist, maxAV);
    if((*chi2vals)[i] <= BADCHI2 && i > 0) {
      while((*chi2vals)[i] <= BADCHI2 && fabs(p[i][i-1] - p[0][i-1]) > 0) {
	p[i][i-1] = p[0][i-1] + (p[i][i-1] - p[0][i-1])/2.;
	(*chi2vals)[i] = GetLnLikelihood(data, p[i], m, isAV, isdist, maxAV);
      }
    }
  }
}

#define DEFAULT_INITIAL_PVECSIZE 10
void incrementparameters_foramoeba(int *Nparameters, int *Ntovary, double ***p, int **ia, int varyparam, double initialval, double stepval)
/* This function sets up the initial simplex for amoeba.
   This is called once for each parameter that one wishes to include in the
   simplex.
   Nparameters - the number of parameters in the simplex before adding the new
                 parameter. On output this will be the number of parameters
		 including the new parameter. Nparameters should be zero on
		 input for the first parameter added to the simplex.
   Ntovary - the total number of parameters that are free to vary. This is
             incremented by the function if varyparam == 1.
   p - pointer to a 2-d array storing the simplex. The function will allocate
       memory for p, and update it, for each parameter added. This array will
       then be passed to amoeba.
   ia - pointer to an integer array storing flags that will indicate whether or
        not the parameter is to be varied. The function will allocate memory for
        the ia array, and update it, for each parameter added. This array will
        be passed to amoeba.
   varyparam - a flag indicating whether or not the new parameter is to be
        varied. 1 - yes, 0 - no.
   initialval - initial value to adopt for the parameter.
   stepval - initial step-size to use for this parameter in defining the
             simplex.
*/
{
  int i,j,k;
  int N;
  N = *Nparameters;
  if(N % DEFAULT_INITIAL_PVECSIZE == 0)
    {
      if(!N)
	{
	  *p = (double **) malloc((DEFAULT_INITIAL_PVECSIZE + 1) * sizeof(double *));
	  *ia = (int *) malloc((DEFAULT_INITIAL_PVECSIZE) * sizeof(int));
	  for(i=0;i< (DEFAULT_INITIAL_PVECSIZE + 1); i++)
	    (*p)[i] = (double *) malloc((DEFAULT_INITIAL_PVECSIZE) * sizeof(double));
	}
      else
	{
	  *p = (double **) realloc((*p), (N + DEFAULT_INITIAL_PVECSIZE + 1) * sizeof(double *));
	  *ia = (int *) realloc((*ia), (N + DEFAULT_INITIAL_PVECSIZE) * sizeof(int));
	  for(i=0;i< (N + 1); i++)
	    (*p)[i] = (double *) realloc((*p)[i], (N + DEFAULT_INITIAL_PVECSIZE) * sizeof(double));
	  for(i=(N+1); i < (N + DEFAULT_INITIAL_PVECSIZE + 1); i++)
	    (*p)[i] = (double *) malloc((N + DEFAULT_INITIAL_PVECSIZE) * sizeof(double));
	}
    }
  if(varyparam)
    (*ia)[N] = 1;
  else
    (*ia)[N] = 0;
  for(i=0;i<(*Ntovary)+1;i++) {
    (*p)[i][N] = initialval;
  }
  if(varyparam) {
    k = (*Ntovary)+1;
    for(j=0;j<N;j++)
      (*p)[k][j] = (*p)[0][j];
    (*p)[k][N] = initialval + stepval;
    (*Ntovary)++;
  }
  (*Nparameters)++;
}

double GetAmoeba_Uncertainties(double *param, int paramindx, int Nparam, double chi20, double errval, _ObsData *data, _ModelData *m, double isAV, double isdist, double maxAV) {
  int i, j, Ntest;
  double chi2test;
  double delp;
  double *ptrial;
  if((ptrial = (double *) malloc(Nparam * sizeof(double))) == NULL) {
    fprintf(stderr,"Memory Allocation Error\n");
    exit(1);
  }
  for(i=0; i < Nparam; i++)
    ptrial[i] = param[i];
  if(!errval) {
    if(param[paramindx] != 0) {
      delp = 0.001*param[paramindx];
    } else
      delp = 0.001;
  }
  else
    delp = errval;
  ptrial[paramindx] = param[paramindx] + delp;
  Ntest = 0;
  chi2test = GetLnLikelihood(data, ptrial, m, isAV, isdist, maxAV);
  if(chi2test == chi20) {
    free(ptrial);
    return 0.0;
  }
  while(Ntest < 100 ? ((chi2test <= BADCHI2) ? 1 : (fabs(chi2test - chi20) < 0.99 ||
					      fabs(chi2test - chi20) > 1.01))
	: 0) {
    while(Ntest < 100 && (chi2test <= BADCHI2)) {
      delp = -delp;
      ptrial[paramindx] = param[paramindx] + delp;
      chi2test = GetLnLikelihood(data, ptrial, m, isAV, isdist, maxAV);
      Ntest++;
      if(!(chi2test <= BADCHI2)) break;
      delp = delp/10.;
      ptrial[paramindx] = param[paramindx] + delp;
      chi2test = GetLnLikelihood(data, ptrial, m, isAV, isdist, maxAV);
      Ntest++;
    }
    if(chi2test == chi20) {
      free(ptrial);
      return 0.0;
    }
    delp = delp/(sqrt(fabs(chi2test - chi20)));
    ptrial[paramindx] = param[paramindx] + delp;
    chi2test = GetLnLikelihood(data, ptrial, m, isAV, isdist, maxAV);
    Ntest++;
  }
  free(ptrial);
  return delp;
}


void amoeba_cleanup(int *Nparameters, int *Ntovary, double ***p, int **ia, double **chi2vals)
{
  int i, k;
  int N, Ns;
  N = *Nparameters;
  if(N <= 0)
    return;
  Ns = ((N - 1) / DEFAULT_INITIAL_PVECSIZE);
  Ns *= DEFAULT_INITIAL_PVECSIZE;
  if(*p != NULL) {
    k = DEFAULT_INITIAL_PVECSIZE + Ns + 1;
    for(i=0; i < k; i++) {
      if((*p)[i] != NULL) {
	free((*p)[i]);
      }
    }
    free((*p));
  }
  if(ia != NULL) {
    if(*ia != NULL)
      free((*ia));
  }
  if(chi2vals != NULL) {
    if(*chi2vals != NULL)
      free((*chi2vals));
  }
}

#undef TINY
#undef SWAP
#undef VARTOOLS_AMOEBA_DEFAULT_MAXSTEPS
#undef DEFAULT_INITIAL_PVECSIZE

/* Functions to generate a guassian random number */
double ran1(void)
{
  return (double) (rand() / (double) RAND_MAX);
}

double gaussian(void)
{
  static int iset = 0;
  static double gset;
  double fac, rsq, v1, v2;

  if(iset == 0) {
    do {
      v1 = 2.0*drand48() - 1.0;
      v2 = 2.0*drand48() - 1.0;
      rsq = v1*v1 + v2*v2;
    } while (rsq >= 1.0 || rsq == 0.0);
    fac = sqrt(-2.0*log(rsq)/rsq);
    gset = v1*fac;
    iset = 1;
    return v2*fac;
  } else {
    iset = 0;
    return gset;
  }
}


/* Functions for performing singular value decomposition and 
   backsubstitution. Taken from Numerical Recipes in C. */
#define ABS_(val) ((val) > 0 ? (val) : -(val))

#define SIGN_(A,B) ((B) >= 0.0 ? ABS_(A) : (-(ABS_(A))))
#define MAX_(A,B) ((A) > (B) ? (A) : (B))
#define MIN_(A,B) ((A) < (B) ? (A) : (B))
#define TINY 1.0e-20
#define TOL 1.0e-9

double pythag(double a, double b)
{
  double absa,absb;
  absa = ABS_(a);
  absb = ABS_(b);
  if(absa > absb) return absa*sqrt(1.0+((absb*absb)/(absa*absa)));
  else return (absb == 0.0 ? 0.0 : absb*sqrt(1.0 + ((absa*absa)/(absb*absb))));
}

void svdcmp(double **a, int m, int n, double *w, double **v)
{
  int flag, i, its, j, jj, k, l, nm;
  double anorm, c, f, g, h, s, scale, x, y, z, *rv1;

  rv1 = (double *) malloc(n*sizeof(double));
  g=scale=anorm=0.0;
  for(i=0;i<n;i++) {
    l=i+1;
    rv1[i]=scale*g;
    g=s=scale=0.0;
    if (i < m) {
      for(k=i;k<m;k++) scale += ABS_(a[k][i]);
      if(scale) {
	for(k=i;k<m;k++) {
	  a[k][i] /= scale;
	  s += a[k][i]*a[k][i];
	}
	f=a[i][i];
	g = -SIGN_(sqrt(s),f);
	h=f*g-s;
	a[i][i]=f-g;
	for (j=l;j<n;j++) {
	  for(s=0.0,k=i;k<m;k++) s += a[k][i]*a[k][j];
	  f=s/h;
	  for (k=i;k<m;k++) a[k][j] += f*a[k][i];
	}
	for (k=i;k<m;k++) a[k][i] *= scale;
      }
    }
    w[i] = scale * g;
    g = s = scale = 0.0;
    if (i < m && i != n - 1) {
      for(k=l;k<n;k++) scale += ABS_(a[i][k]);
      if(scale) {
	for(k=l;k<n;k++) {
	  a[i][k] /= scale;
	  s += a[i][k]*a[i][k];
	}
	f = a[i][l];
	g = -SIGN_(sqrt(s),f);
	h = f*g-s;
	a[i][l] = f-g;
	for (k=l;k<n;k++) rv1[k] = a[i][k]/h;
	for (j=l;j<m;j++) {
	  for (s=0.0,k=l;k<n;k++) s += a[j][k]*a[i][k];
	  for (k=l;k<n;k++) a[j][k] += s*rv1[k];
	}
	for(k=l;k<n;k++) a[i][k] *= scale;
      }
    }
    anorm=MAX_(anorm,(ABS_(w[i])+ABS_(rv1[i])));
  }
  for(i=n-1;i>=0;i--) {
    if (i < n-1) {
      if(g) {
	for (j=l;j<n;j++)
	  v[j][i]=(a[i][j]/a[i][l])/g;
	for(j=l;j<n;j++) {
	  for(s=0.0,k=l;k<n;k++) s += a[i][k]*v[k][j];
	  for(k=l;k<n;k++) v[k][j] += s*v[k][i];
	}
      }
      for(j=l;j<n;j++) v[i][j] = v[j][i] = 0.0;
    }
    v[i][i] = 1.0;
    g = rv1[i];
    l = i;
  }
  for(i=MIN_(m,n)-1;i>=0;i--) {
    l = i+1;
    g = w[i];
    for (j=l;j < n; j++) a[i][j] = 0.0;
    if (g) {
      g = 1.0/g;
      for (j=l;j<n;j++) {
	for (s=0.0,k=l;k<m;k++) s += a[k][i] * a[k][j];
	f = (s/a[i][i]) * g;
	for (k=i;k<m;k++) a[k][j] += f*a[k][i];
      }
      for (j=i;j<m;j++) a[j][i] *= g;
    } else for (j=i;j<m;j++) a[j][i] = 0.0;
    ++a[i][i];
  }
  for (k=n-1;k>=0;k--) {
    for (its=1;its<=30;its++) {
      flag = 1;
      for (l=k;l>=0;l--) {
	nm=l-1;
	if ((double)(ABS_(rv1[l])+anorm) == anorm) {
	  flag = 0;
	  break;
	}
	if ((double)(ABS_(w[nm])+anorm) == anorm) break;
      }
      if(flag) {
	c=0.0;
	s=1.0;
	for (i=l;i<=k;i++) {
	  f = s*rv1[i];
	  rv1[i] = c*rv1[i];
	  if ((double) (ABS_(f)+anorm) == anorm) break;
	  g = w[i];
	  h = pythag(f,g);
	  w[i] = h;
	  h = 1.0/h;
	  c = g*h;
	  s = -f*h;
	  for (j=0;j<m;j++) {
	    y = a[j][nm];
	    z = a[j][i];
	    a[j][nm] = y*c+z*s;
	    a[j][i] = z*c-y*s;
	  }
	}
      }
      z = w[k];
      if (l == k) {
	if (z < 0.0) {
	  w[k] = -z;
	  for (j=0;j<n;j++) v[j][k] = -v[j][k];
	}
	break;
      }
      if (its == 30)
	{
	  fprintf(stderr,"Error: No convergence in 30 svdcmp iterations\n");
	  exit(5);
	}
      x = w[l];
      nm = k-1;
      y=w[nm];
      g=rv1[nm];
      h=rv1[k];
      f = ((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
      g = pythag(f,1.0);
      f = ((x-z)*(x+z)+h*((y/(f+SIGN_(g,f)))-h))/x;
      c = s = 1.0;
      for (j=l;j<=nm;j++) {
	i = j+1;
	g = rv1[i];
	y=w[i];
	h = s*g;
	g = c*g;
	z = pythag(f,h);
	rv1[j] = z;
	c = f/z;
	s = h/z;
	f = x*c+g*s;
	g = g*c-x*s;
	h=y*s;
	y *= c;
	for (jj=0;jj<n;jj++) {
	  x = v[jj][j];
	  z = v[jj][i];
	  v[jj][j] = x*c + z*s;
	  v[jj][i] = z*c - x*s;
	}
	z = pythag(f,h);
	w[j] = z;
	if (z) {
	  z = 1.0/z;
	  c = f*z;
	  s = h*z;
	}
	f = c*g + s*y;
	x = c*y - s*g;
	for (jj=0;jj<m;jj++) {
	  y = a[jj][j];
	  z = a[jj][i];
	  a[jj][j] = y*c + z*s;
	  a[jj][i] = z*c - y*s;
	}
      }
      rv1[l] = 0.0;
      rv1[k] = f;
      w[k] = x;
    }
  }
  free(rv1);
}


void svbksb(double **u, double *w, double **v, int m, int n, double *b, double *x)
{
  int jj,j,i;
  double s,*tmp;
  tmp = (double *) malloc(n*sizeof(double));
  for(j=0;j<n;j++) {
    s = 0.0;
    if(w[j]) {
      for(i=0;i<m;i++) {if(!isnan(b[i])) s += u[i][j]*b[i];}
      s /= w[j];
    }
    tmp[j] = s;
  }
  for(j=0;j<n;j++) {
    s = 0.0;
    for(jj=0;jj<n;jj++) s += v[j][jj]*tmp[jj];
    x[j]=s;
  }
  free(tmp);
}

#undef ABS_

#undef SIGN_
#undef MAX_
#undef MIN_
#undef TINY 
#undef TOL 

void InitializeEmpiricalModelParams(_EmpiricalModelParams *mod, int isrand) {
  double r;
  int i;
  do {
    for(i=0; i < NEMPIRICALMODELPARAMS; i++)
      mod->v[i] = ModelMedian[i];
    if(isrand) {
      for(i=0; i < NEMPIRICALMODELPARAMS; i++) {
	r = gaussian();
	mod->v[i] += r*ModelStddev[i];
      }
    }
  } while(mod->v[INDEX_Sr] < TINYERRTERM || mod->v[INDEX_Sm] < TINYERRTERM );
}

void EvaluateEmpiricalModel(_ModelData *m, _EmpiricalModelParams *mod, double logrho, double dlR, double dMbol) {
  double logrho2, logrho3, lM, lM2, lM3, lT, lT2;
  double Crho = 0.1496766;
  double Cmbol = 42.227;

  m->logrho = logrho;
  m->dlR = dlR;
  m->dMbol = dMbol;

  logrho2 = logrho*logrho; logrho3 = logrho2*logrho;
  m->lR = 
    mod->v[INDEX_a0]*(0.131674 - 0.605624*logrho + 0.739127*logrho2 
		      - 0.263765*logrho3)+
    mod->v[INDEX_a1]*(-0.494952 + 0.614479*logrho + 0.437886*logrho2 
		      -0.430922*logrho3)+
    mod->v[INDEX_a2]*(-0.765425 - 0.291096*logrho + 0.099548*logrho2 
		      +0.565224*logrho3)+
    mod->v[INDEX_a3]*(0.389626 + 0.413398*logrho + 0.502032*logrho2 
		      +0.652117*logrho3)+
    dlR;
  m->lM = logrho - Crho + 3.0*m->lR;
  lM = m->lM; lM2 = lM*lM; lM3 = lM2*lM;
  m->Mbol =
    mod->v[INDEX_b0]*(-0.042114 - 0.349391*lM - 0.787638*lM2 - 0.505746*lM3) + 
    mod->v[INDEX_b1]*(0.170316 + 0.723492*lM + 0.095153*lM2 - 0.662192*lM3)+
    mod->v[INDEX_b2]*(-0.309788 - 0.505183*lM + 0.591641*lM2 - 0.546610*lM3)+
    mod->v[INDEX_b3]*(-0.934479 + 0.315080*lM - 0.143296*lM2 + 0.083309*lM3)+
    dMbol;
  m->lTeff = (Cmbol - 5*m->lR - m->Mbol)/10.0;
  lT = m->lTeff; lT2 = lT*lT;
  m->BCV = 
    mod->v[INDEX_aBCV0]*(0.871464 - 0.485768*lT + 0.067675*lT2) + 
    mod->v[INDEX_aBCV1]*(0.484681 + 0.831841*lT - 0.270418*lT2) + 
    mod->v[INDEX_aBCV2]*(0.075065 + 0.268461*lT + 0.960361*lT2);
  m->BCJ = 
    mod->v[INDEX_aBCJ0]*(0.871464 - 0.485768*lT + 0.067675*lT2) + 
    mod->v[INDEX_aBCJ1]*(0.484681 + 0.831841*lT - 0.270418*lT2) + 
    mod->v[INDEX_aBCJ2]*(0.075065 + 0.268461*lT + 0.960361*lT2);
  m->BCH = 
    mod->v[INDEX_aBCH0]*(0.871464 - 0.485768*lT + 0.067675*lT2) + 
    mod->v[INDEX_aBCH1]*(0.484681 + 0.831841*lT - 0.270418*lT2) + 
    mod->v[INDEX_aBCH2]*(0.075065 + 0.268461*lT + 0.960361*lT2);
  m->BCK = 
    mod->v[INDEX_aBCK0]*(0.871464 - 0.485768*lT + 0.067675*lT2) + 
    mod->v[INDEX_aBCK1]*(0.484681 + 0.831841*lT - 0.270418*lT2) + 
    mod->v[INDEX_aBCK2]*(0.075065 + 0.268461*lT + 0.960361*lT2);
  m->absVmag = m->Mbol + m->BCV;
  m->absJmag = m->Mbol + m->BCJ;
  m->absHmag = m->Mbol + m->BCH;
  m->absKmag = m->Mbol + m->BCK;  
}

/* Adds uncertainty in the empirical model parameters to the likelihood 
   function */
double EmpiricalModelLikelihood(_EmpiricalModelParams *mod) {
  double L;
  double a;
  int j;
  a = (mod->v[0] - ModelMedian[0])/ModelStddev[0];
  L = a*a;
  for(j = 1; j < NEMPIRICALMODELPARAMS; j++) {
    a = (mod->v[j] - ModelMedian[j])/ModelStddev[j];
    L += a*a;
  }
  return L;
}


/* Takes as input parameters:
   Input Observation struction,
   lrho, dlR, dMbol, AV, distmod, and the set of empirical model parameters
   Outputs:
   -2*ln(likelihood)
*/
double GetLnLikelihood(_ObsData *data, double *paramvec, _ModelData *m, double isAV, double isdist, double maxAV) {
  double L, mval, diff;
  double lrho, dlR, dMbol;
  _EmpiricalModelParams *mod;
  double distmod; double AV;
  int i, j;
  lrho = paramvec[0];
  dlR = paramvec[1];
  dMbol = paramvec[2];
  mod = (_EmpiricalModelParams *) (&(paramvec[3]));
  j = NEMPIRICALMODELPARAMS + 3;
  if(isdist) {distmod = paramvec[j]; j++;} else distmod=0.;
  if(isAV) {AV = paramvec[j]; j++;} else AV = 0.;
  if(AV < 0) return BADCHI2;
  if(maxAV >= 0 && AV > maxAV) return BADCHI2;
  L = 0.0;
  if(lrho < MIN_LRHO_LIMIT || lrho > MAX_LRHO_LIMIT ) return BADCHI2;
  EvaluateEmpiricalModel(m, mod, lrho, dlR, dMbol);
  m->AV = AV; m->distmod = distmod;
  for(i=0; i < data->Nobs; i++) {
    switch(data->obstype[i]) {
    case OBS_TYPE_MASS:
      mval = exp(M_LN10*m->lM);
      break;
    case OBS_TYPE_RADIUS:
      mval = exp(M_LN10*m->lR);
      break;
    case OBS_TYPE_DENSITY:
      mval = exp(M_LN10*lrho);
      break;
    case OBS_TYPE_LOGDENSITY:
      mval = lrho;
      break;
    case OBS_TYPE_TEFF:
      mval = exp(M_LN10*m->lTeff);
      break;
    case OBS_TYPE_VMAG:
      mval = m->absVmag + distmod + AV;
      break;
    case OBS_TYPE_JMAG:
      mval = m->absJmag + distmod + AV*REDLAW_AJ;
      break;
    case OBS_TYPE_HMAG:
      mval = m->absHmag + distmod + AV*REDLAW_AH;
      break;
    case OBS_TYPE_KMAG:
      mval = m->absKmag + distmod + AV*REDLAW_AK;
      break;
    case OBS_TYPE_ABSVMAG:
      mval = m->absVmag;
      break;
    case OBS_TYPE_ABSJMAG:
      mval = m->absJmag;
      break;
    case OBS_TYPE_ABSHMAG:
      mval = m->absHmag;
      break;
    case OBS_TYPE_ABSKMAG:
      mval = m->absKmag;
      break;
    case OBS_TYPE_VJCOLOR:
      mval = m->absVmag - m->absJmag + AV*(1 - REDLAW_AJ);
      break;
    case OBS_TYPE_VHCOLOR:
      mval = m->absVmag - m->absHmag + AV*(1 - REDLAW_AH);
      break;
    case OBS_TYPE_VKCOLOR:
      mval = m->absVmag - m->absKmag + AV*(1 - REDLAW_AK);
      break;
    case OBS_TYPE_JHCOLOR:
      mval = m->absJmag - m->absHmag + AV*(REDLAW_AJ - REDLAW_AH);
      break;
    case OBS_TYPE_JKCOLOR:
      mval = m->absJmag - m->absKmag + AV*(REDLAW_AJ - REDLAW_AK);
      break;
    case OBS_TYPE_HKCOLOR:
      mval = m->absHmag - m->absKmag + AV*(REDLAW_AH - REDLAW_AK);
      break;
    case OBS_TYPE_DISTMOD:
      mval = distmod;
      break;
    case OBS_TYPE_DIST:
      mval = exp(M_LN10*((distmod + 5.0)/5.0));
      break;
    case OBS_TYPE_PARALLAX:
      mval = 1000.0/exp(M_LN10*((distmod + 5.0)/5.0));
      break;
    case OBS_TYPE_ANGDIAM:
      mval = 0.009298438*exp(M_LN10*m->lR)*1000.0/exp(M_LN10*((distmod + 5.0)/5.0));
      break;
    case OBS_TYPE_AV:
      mval = AV;
      break;
    case OBS_TYPE_EBV:
      mval = AV*(REDLAW_AB - 1.0);
      break;
    case OBS_TYPE_LOGG:
      mval = LOGGSUN + m->lM - 2.0*m->lR;
      break;
    case OBS_TYPE_LOGGSOLAR:
      mval = m->lM - 2.0*m->lR;
      break;
    case OBS_TYPE_DENSITYSOLAR:
      mval = exp(M_LN10*(m->lM - 3.0*m->lR));
      break;
    case OBS_TYPE_LOGDENSITYSOLAR:
      mval = m->lM - 3.0*m->lR;
      break;
    default:
      return BADCHI2;
      break;
    }
    if(data->islower[i]) {
      if(mval < data->lowerlim[i]) {
	return BADCHI2;
      }
    }
    if(data->isupper[i]) {
      if(mval > data->upperlim[i]) {
	return BADCHI2;
      }
    }
    diff = (mval - data->obsval[i])/data->errval[i];
    L += (diff*diff);
  }
  L += (dlR*dlR/mod->v[INDEX_Sr]/mod->v[INDEX_Sr] + 
	log(mod->v[INDEX_Sr]*mod->v[INDEX_Sr]));
  L += (dMbol*dMbol/mod->v[INDEX_Sm]/mod->v[INDEX_Sm] + 
	log(mod->v[INDEX_Sm]*mod->v[INDEX_Sm]));
  L += EmpiricalModelLikelihood(mod);
  return L;
}

void parseobsdata(char *argv0, char *argvpar, _ObsData *obs, char *argv, int obstype) {
  int i, j, k;
  char *varstrcpy;
  int l;
  if(obs->Nobs == 0) {
    if((obs->obsval = (double *) malloc(sizeof(double))) == NULL ||
       (obs->errval = (double *) malloc(sizeof(double))) == NULL ||
       (obs->obstype = (int *) malloc(sizeof(int))) == NULL ||
       (obs->isupper = (char *) malloc(1)) == NULL ||
       (obs->islower = (char *) malloc(1)) == NULL ||
       (obs->upperlim = (double *) malloc(sizeof(double))) == NULL ||
       (obs->lowerlim = (double *) malloc(sizeof(double))) == NULL) {
      fprintf(stderr,"Memory Allocation Error\n");
      exit(1);
    }
  }
  else {
    if((obs->obsval = (double *) realloc(obs->obsval, (obs->Nobs + 1)*sizeof(double))) == NULL ||
       (obs->errval = (double *) realloc(obs->errval, (obs->Nobs + 1)*sizeof(double))) == NULL ||
       (obs->obstype = (int *) realloc(obs->obstype, (obs->Nobs + 1)*sizeof(int))) == NULL ||
       (obs->isupper = (char *) realloc(obs->isupper, (obs->Nobs + 1))) == NULL ||
       (obs->islower = (char *) realloc(obs->islower, (obs->Nobs + 1))) == NULL ||
       (obs->upperlim = (double *) realloc(obs->upperlim, (obs->Nobs + 1)*sizeof(double))) == NULL ||
       (obs->lowerlim = (double *) realloc(obs->lowerlim, (obs->Nobs + 1)*sizeof(double))) == NULL) {
      fprintf(stderr,"Memory Allocation Error\n");
      exit(1);
    }
  }
  
  obs->obstype[obs->Nobs] = obstype;
  obs->islower[obs->Nobs] = 0;
  obs->isupper[obs->Nobs] = 0;

  l = strlen(argv)+1;
  varstrcpy = malloc(l);
  sprintf(varstrcpy,"%s",argv);

  i = 0;
  while(varstrcpy[i] != '\0' && varstrcpy[i] != ':' && varstrcpy[i] != '[') {
    i++;
  }
  if(varstrcpy[i] != ':') {
    fprintf(stderr,"Error parsing %s %s\nThe expected format is %s obsval:errval[lowerlim:upperlim] where [lowerlim:upperlim] is optional.\n", argvpar, argv, argvpar);
    exit(1);
  }
  varstrcpy[i] = '\0';
  j = sscanf(varstrcpy,"%lf",&(obs->obsval[obs->Nobs]));
  if(!j) {
    fprintf(stderr,"Error parsing %s %s\nThe expected format is %s obsval:errval[lowerlim:upperlim] where [lowerlim:upperlim] is optional.\n", argvpar, argv, argvpar);
    exit(1);
  }

  k = i+1;
  i = i+1;
  while(varstrcpy[i] != '\0' && varstrcpy[i] != ':' && varstrcpy[i] != '[') {
    i++;
  }
  if(varstrcpy[i] != '[') {
    varstrcpy[i] = '\0';
    j = sscanf(&(varstrcpy[k]),"%lf",&(obs->errval[obs->Nobs]));
    obs->Nobs = obs->Nobs + 1;
    return;
  }

  varstrcpy[i] = '\0';
  j = sscanf(&(varstrcpy[k]),"%lf",&(obs->errval[obs->Nobs]));

  k = i+1;
  i = i+1;
  while(varstrcpy[i] != '\0' && varstrcpy[i] != ':' && varstrcpy[i] != ']') {
    i++;
  }
  if(varstrcpy[i] != ':') {
    varstrcpy[i] = '\0';
    if(i != k) {
      j = sscanf(&(varstrcpy[k]),"%lf",&(obs->lowerlim[obs->Nobs]));
      if(j)
	obs->islower[obs->Nobs] = 1;
    }
    obs->Nobs = obs->Nobs + 1;
    return;
  }
  
  varstrcpy[i] = '\0';
  if(i != k) {
    j = sscanf(&(varstrcpy[k]),"%lf",&(obs->lowerlim[obs->Nobs]));
    if(j)
      obs->islower[obs->Nobs] = 1;
  }

  k = i+1;
  i = i+1;
  while(varstrcpy[i] != '\0' && varstrcpy[i] != ']') {
    i++;
  }

  varstrcpy[i] = '\0';
  if(i != k) {
    j = sscanf(&(varstrcpy[k]),"%lf",&(obs->upperlim[obs->Nobs]));
    if(j)
      obs->isupper[obs->Nobs] = 1;
  }

  obs->Nobs += 1;

}

/* Fit a polynomial between x and y data, evaluate it at a specified value
   of x, and return the result */
#define SVDTOLERANCE 1.0e-9
double FitAndEvaluatePoly(double *x, double *y, int N, double xeval, int order)
{
  int i, j;
  int ncomp;
  double yeval;
  double **Amatrix;
  double *Bvector;
  double *Avector;
  double **v, *w, *wti;
  double term, wmax, thresh;
  ncomp = order + 1;
  if((Amatrix = (double **) malloc(N*sizeof(double *))) == NULL ||
     (Bvector = (double *) malloc(N*sizeof(double))) == NULL ||
     (Avector = (double *) malloc(ncomp*sizeof(double))) == NULL ||
     (v = (double **) malloc(ncomp*sizeof(double *))) == NULL ||
     (w = (double *) malloc(ncomp*sizeof(double))) == NULL ||
     (wti = (double *) malloc(ncomp*sizeof(double))) == NULL) {
    fprintf(stderr,"Memory Allocation Error\n");
    exit(1);
  }
  for(i=0; i < N; i++) {
    if((Amatrix[i] = (double *) malloc(ncomp * sizeof(double))) == NULL)
      {
	fprintf(stderr,"Memory Allocation Error\n");
	exit(1);
      }
  }
  for(i=0; i < ncomp; i++) {
    if((v[i] = (double *) malloc(ncomp * sizeof(double))) == NULL) {
      fprintf(stderr,"Memory Allocation Error\n");
      exit(1);
    }
  }
  /* Fill out the design matrix */
  for(i=0; i < N; i++) {
    Bvector[i] = y[i];
    Amatrix[i][0] = 1.;
    term = 1.;
    for(j=1; j <= order; j++) {
      term *= x[i];
      Amatrix[i][j] = term;
    }
  }
  /* Run singular value decomposition to fit the model */
  svdcmp(Amatrix,N,ncomp,w,v);
  wmax = 0.;
  for(j=0; j < ncomp; j++) {
    if(w[j] > wmax) wmax = w[j];
  }
  thresh = SVDTOLERANCE * wmax;
  for(j=0; j < ncomp; j++) {
    if(w[j] < thresh) w[j] = 0.0;
  }
  svbksb(Amatrix,w,v,N,ncomp,Bvector,Avector);

  /* Evaluate the polynomial */
  yeval = Avector[0];
  term = 1;
  for(j=1; j <= order; j++) {
    term *= xeval;
    yeval += term*Avector[j];
  }
  
  /* Clean up the allocated memory and return the result */
  for(i=0; i < ncomp; i++) free(v[i]);
  for(i=0; i < N; i++) free(Amatrix[i]);
  free(Amatrix);
  free(Bvector);
  free(Avector);
  free(v);
  free(w);
  free(wti);

  return yeval;
}

/* This function determines the initial density estimate from a value 
 it does this by evaluating the empirical model over a range of densities
 and then fitting a relation between density and the value.  While this is
 slower than hard-coding in inverse relations, this should make the code
 easier to update when the empirical model is revised.
*/

#define NSTDDEV_SIM 100
#define NRHOSTEP 40

void GetDensityFromValue(double *logrhoout, double *logrhoerrout, int obstype, double valin, double errin) {
  void PrintDEMCChain(FILE *outfile, _ModelData *m, double *l, int N);
  double logrhostep[NRHOSTEP] = {0.000000, 0.050000, 0.100000, 0.150000, 0.200000, 0.250000, 0.300000, 0.350000, 0.400000, 0.450000, 0.500000, 0.550000, 0.600000, 0.650000, 0.700000, 0.750000, 0.800000, 0.850000, 0.900000, 0.950000, 1.000000, 1.050000, 1.100000, 1.150000, 1.200000, 1.250000, 1.300000, 1.350000, 1.400000, 1.450000, 1.500000, 1.550000, 1.600000, 1.650000, 1.700000, 1.750000, 1.800000, 1.850000, 1.900000, 1.950000};
  int Nrhostep = NRHOSTEP;
  _ModelData m[NRHOSTEP];
  _EmpiricalModelParams mod;
  double xvals[NRHOSTEP];
  double vals_stddev[NSTDDEV_SIM];
  int i, j;
  int Nfitstart, Nfitstop;

  /* Set the empirical model parameters to the median values */
  InitializeEmpiricalModelParams(&mod, 0);

  /* Evaluate the empirical model over a range of densities */
  for(i=0; i < Nrhostep; i++) {
    EvaluateEmpiricalModel(&(m[i]), &mod, logrhostep[i], 0.0, 0.0);
  }
  
  /* Extract the parameter that we are using for estimating the density
     from each model evaluation */
  for(i=0; i < Nrhostep; i++) {
    switch(obstype) {
    case OBS_TYPE_RADIUS:
      Nfitstart = 0; Nfitstop = Nrhostep-1;
      xvals[i] = exp(M_LN10*m[i].lR);
      break;
    case OBS_TYPE_MASS:
      Nfitstart = 2; Nfitstop = Nrhostep-1;
      xvals[i] = exp(M_LN10*m[i].lM);
      break;
    case OBS_TYPE_TEFF:
      Nfitstart = 4; Nfitstop = Nrhostep-1;
      xvals[i] = exp(M_LN10*m[i].lTeff);
      break;
    case OBS_TYPE_ABSVMAG:
      Nfitstart = 3; Nfitstop = Nrhostep-1;
      xvals[i] = m[i].absVmag;
      break;
    case OBS_TYPE_ABSJMAG:
      Nfitstart = 0; Nfitstop = Nrhostep-1;
      xvals[i] = m[i].absJmag;
      break;
    case OBS_TYPE_ABSHMAG:
      Nfitstart = 1; Nfitstop = Nrhostep-1;
      xvals[i] = m[i].absHmag;
      break;
    case OBS_TYPE_ABSKMAG:
      Nfitstart = 1; Nfitstop = Nrhostep-1;
      xvals[i] = m[i].absKmag;
      break;
    case OBS_TYPE_VJCOLOR:
      Nfitstart = 4; Nfitstop = Nrhostep-1;
      xvals[i] = m[i].absVmag - m[i].absJmag;
      break;
    case OBS_TYPE_VHCOLOR:
      Nfitstart = 4; Nfitstop = Nrhostep-1;
      xvals[i] = m[i].absVmag - m[i].absHmag;
      break;
    case OBS_TYPE_VKCOLOR:
      Nfitstart = 4; Nfitstop = Nrhostep-1;
      xvals[i] = m[i].absVmag - m[i].absKmag;
      break;
    case OBS_TYPE_JHCOLOR:
      Nfitstart = 4; Nfitstop = Nrhostep-17;
      xvals[i] = m[i].absJmag - m[i].absHmag;
      break;
    case OBS_TYPE_JKCOLOR:
      if(m[i].lTeff > 3.6) {
	xvals[i] = -12.1807*m[i].lTeff*m[i].lTeff + 86.3564*m[i].lTeff - 152.202;
      } else {
	xvals[i] = m[i].absJmag - m[i].absKmag;
      }
      Nfitstart = 4; Nfitstop = Nrhostep-1;
      break;
    case OBS_TYPE_HKCOLOR:
      Nfitstart = 4; Nfitstop = Nrhostep-1;
      xvals[i] = m[i].absHmag - m[i].absKmag;
      break;
    case OBS_TYPE_LOGG:
      Nfitstart = 0; Nfitstop = Nrhostep-1;
      xvals[i] = LOGGSUN + m[i].lM - 2.0*m[i].lR;
      break;
    case OBS_TYPE_LOGGSOLAR:
      Nfitstart = 0; Nfitstop = Nrhostep-1;
      xvals[i] = m[i].lM - 2.0*m[i].lR;
      break;
    default:
      fprintf(stderr,"Error: Unsupported initialization parameter. There is a bug in the code. Contact the developer giving the command line used, and this message.\n");
      exit(1);
    }
  }
  
  /* Fit a polynomial between the estimator and the density, and then
     evaluate it at the observed value for the estimator */
  *logrhoout = FitAndEvaluatePoly(&(xvals[Nfitstart]), &(logrhostep[Nfitstart]), (Nfitstop - Nfitstart + 1), valin, 5);

  /* To initialize the density uncertainty repeat the above 100 times and 
     find the scatter in the density */
  for(j=0; j < NSTDDEV_SIM; j++) {
    /* Set the empirical model parameters, this time allowing them to
     vary*/
    InitializeEmpiricalModelParams(&mod, 1);

    /* Evaluate the empirical model over a range of densities */
    for(i=0; i < Nrhostep; i++) {
      EvaluateEmpiricalModel(&(m[i]), &mod, logrhostep[i], (mod.v[INDEX_Sr]*gaussian()), (mod.v[INDEX_Sm]*gaussian()));
    }
  
    /* Extract the parameter that we are using for estimating the density
       from each model evaluation */
    for(i=0; i < Nrhostep; i++) {
      switch(obstype) {
      case OBS_TYPE_RADIUS:
	Nfitstart = 0; Nfitstop = Nrhostep-1;
	xvals[i] = exp(M_LN10*m[i].lR);
	break;
      case OBS_TYPE_MASS:
	Nfitstart = 2; Nfitstop = Nrhostep-1;
	xvals[i] = exp(M_LN10*m[i].lM);
	break;
      case OBS_TYPE_TEFF:
	Nfitstart = 4; Nfitstop = Nrhostep-1;
	xvals[i] = exp(M_LN10*m[i].lTeff);
	break;
      case OBS_TYPE_ABSVMAG:
	Nfitstart = 3; Nfitstop = Nrhostep-1;
	xvals[i] = m[i].absVmag;
	break;
      case OBS_TYPE_ABSJMAG:
	Nfitstart = 0; Nfitstop = Nrhostep-1;
	xvals[i] = m[i].absJmag;
	break;
      case OBS_TYPE_ABSHMAG:
	Nfitstart = 1; Nfitstop = Nrhostep-1;
	xvals[i] = m[i].absHmag;
	break;
      case OBS_TYPE_ABSKMAG:
	Nfitstart = 1; Nfitstop = Nrhostep-1;
	xvals[i] = m[i].absKmag;
	break;
      case OBS_TYPE_VJCOLOR:
	Nfitstart = 4; Nfitstop = Nrhostep-1;
	xvals[i] = m[i].absVmag - m[i].absJmag;
	break;
      case OBS_TYPE_VHCOLOR:
	Nfitstart = 4; Nfitstop = Nrhostep-1;
	xvals[i] = m[i].absVmag - m[i].absHmag;
	break;
      case OBS_TYPE_VKCOLOR:
	Nfitstart = 4; Nfitstop = Nrhostep-1;
	xvals[i] = m[i].absVmag - m[i].absKmag;
	break;
      case OBS_TYPE_JHCOLOR:
	Nfitstart = 4; Nfitstop = Nrhostep-17;
	xvals[i] = m[i].absJmag - m[i].absHmag;
	break;
      case OBS_TYPE_JKCOLOR:
	Nfitstart = 4; Nfitstop = Nrhostep-1;
	if(m[i].lTeff > 3.6) {
	  xvals[i] = -12.1807*m[i].lTeff*m[i].lTeff + 86.3564*m[i].lTeff - 152.202;
	} else {
	  xvals[i] = m[i].absJmag - m[i].absKmag;
	}
	break;
      case OBS_TYPE_HKCOLOR:
	Nfitstart = 4; Nfitstop = Nrhostep-1;
	xvals[i] = m[i].absHmag - m[i].absKmag;
	break;
      case OBS_TYPE_LOGG:
	Nfitstart = 0; Nfitstop = Nrhostep-1;
	xvals[i] = LOGGSUN + m[i].lM - 2.0*m[i].lR;
	break;
      case OBS_TYPE_LOGGSOLAR:
	Nfitstart = 0; Nfitstop = Nrhostep-1;
	xvals[i] = m[i].lM - 2.0*m[i].lR;
	break;
      default:
	fprintf(stderr,"Error: Unsupported initialization parameter. There is a bug in the code. Contact the developer giving the command line used, and this message.\n");
	exit(1);
      }
    }
    
    /* Fit a polynomial between the estimator and the density, and then
       evaluate it at the observed value for the estimator */
    vals_stddev[j] = FitAndEvaluatePoly(&(xvals[Nfitstart]), &(logrhostep[Nfitstart]), (Nfitstop - Nfitstart + 1), (valin + errin*gaussian()), 5);
    
  }
  *logrhoerrout = ComputeStddev(vals_stddev,NSTDDEV_SIM);
}


/* This function determines the initial estimates for the distance modulus and AV */
void GuessInitialDistAV(double *initdistmod,double *initdistmoderr,double *initAV,double *initAVerr,int isdistmod,int isAV,_ObsData *data,double maxAV,double logrhoinit)
{
  int i;
  int founddist = 0;
  int Nmag = 0, Ncolor = 0;
  _ModelData mdata;
  _EmpiricalModelParams mod;
  double AV, distmod, eAV, edistmod;
  double a1_ = 0., a2_ = 0., b1_ = 0., b2_ = 0., c1_ = 0., c2_ = 0.;
  double Vmag, Jmag, Hmag, Kmag, eVmag, eJmag, eHmag, eKmag;
  double VJ, VH, VK, JH, JK, HK, eVJ, eVH, eVK, eJH, eJK, eHK;
  char isV = 0, isJ = 0, isH = 0, isK = 0, isVJ = 0, isVH = 0, isVK = 0, isJH = 0, isJK = 0, isHK = 0;
  if(isdistmod) {
    for(i=0; i < data->Nobs; i++) {
      if(data->obstype[i] == OBS_TYPE_DIST) {
	*initdistmod = 5*log(data->obsval[i])/M_LN10 - 5.0;
	*initdistmoderr = 5*data->errval[i]/data->obsval[i]/M_LN10;
	distmod = *initdistmod;
	edistmod = *initdistmoderr;
	founddist = 1;
      }
      else if(data->obstype[i] == OBS_TYPE_DISTMOD) {
	*initdistmod = data->obsval[i];
	*initdistmoderr = data->errval[i];
	distmod = *initdistmod;
	edistmod = *initdistmoderr;
	founddist = 1;
      }
      else if(data->obstype[i] == OBS_TYPE_PARALLAX) {
	*initdistmod = 5*log((1000.0/data->obsval[i]))/M_LN10 - 5.0;
	*initdistmoderr = 5*data->errval[i]/data->obsval[i]/M_LN10;
	distmod = *initdistmod;
	edistmod = *initdistmoderr;
	founddist = 1;
      }
      else if(data->obstype[i] == OBS_TYPE_ANGDIAM && !founddist) {
	InitializeEmpiricalModelParams(&mod, 0);
	EvaluateEmpiricalModel(&mdata, &mod, logrhoinit, 0., 0.);

	*initdistmod = 5*log(9.298438*exp(M_LN10*mdata.lR)/data->obsval[i])/M_LN10 - 5.0;
	*initdistmoderr = 5*data->errval[i]/data->obsval[i]/M_LN10;

	distmod = *initdistmod;
	edistmod = *initdistmoderr;
	founddist = 1;
      }
    }
  }
  
  /* Get all the magnitudes and colors from the observations */
  for(i=0; i < data->Nobs; i++) {
    if(data->obstype[i] == OBS_TYPE_VMAG) {
      if(!isV) { 
	Vmag = data->obsval[i]; eVmag = data->errval[i]; isV = 1;
	Nmag++;
      } else if(data->errval[i] < eVmag) {
	Vmag = data->obsval[i]; eVmag = data->errval[i]; isV = 1;
      }
    }
    else if(data->obstype[i] == OBS_TYPE_JMAG) {
      if(!isJ) { 
	Jmag = data->obsval[i]; eJmag = data->errval[i]; isJ = 1;
	Nmag++;
      } else if(data->errval[i] < eJmag) {
	Jmag = data->obsval[i]; eJmag = data->errval[i]; isJ = 1;
      }
    }
    else if(data->obstype[i] == OBS_TYPE_HMAG) {
      if(!isH) { 
	Hmag = data->obsval[i]; eHmag = data->errval[i]; isH = 1;
	Nmag++;
      } else if(data->errval[i] < eHmag) {
	Hmag = data->obsval[i]; eHmag = data->errval[i]; isH = 1;
      }
    }
    else if(data->obstype[i] == OBS_TYPE_KMAG) {
      if(!isK) { 
	Kmag = data->obsval[i]; eKmag = data->errval[i]; isK = 1;
	Nmag++;
      } else if(data->errval[i] < eKmag) {
	Kmag = data->obsval[i]; eKmag = data->errval[i]; isK = 1;
      }
    }
    else if(data->obstype[i] == OBS_TYPE_VJCOLOR) {
      if(!isVJ) { 
	VJ = data->obsval[i]; eVJ = data->errval[i]; isVJ = 1;
	Ncolor++;
      } else if(data->errval[i] < eVJ) {
	VJ = data->obsval[i]; eVJ = data->errval[i]; isVJ = 1;
      }
    }
    else if(data->obstype[i] == OBS_TYPE_VHCOLOR) {
      if(!isVH) { 
	VH = data->obsval[i]; eVH = data->errval[i]; isVH = 1;
	Ncolor++;
      } else if(data->errval[i] < eVH) {
	VH = data->obsval[i]; eVH = data->errval[i]; isVH = 1;
      }
    }
    else if(data->obstype[i] == OBS_TYPE_VKCOLOR) {
      if(!isVK) { 
	VK = data->obsval[i]; eVK = data->errval[i]; isVK = 1;
	Ncolor++;
      } else if(data->errval[i] < eVK) {
	VK = data->obsval[i]; eVK = data->errval[i]; isVK = 1;
      }
    }
    else if(data->obstype[i] == OBS_TYPE_JHCOLOR) {
      if(!isJH) { 
	JH = data->obsval[i]; eJH = data->errval[i]; isJH = 1;
	Ncolor++;
      } else if(data->errval[i] < eJH) {
	JH = data->obsval[i]; eJH = data->errval[i]; isJH = 1;
      }
    }
    else if(data->obstype[i] == OBS_TYPE_JKCOLOR) {
      if(!isJK) { 
	JK = data->obsval[i]; eJK = data->errval[i]; isJK = 1;
	Ncolor++;
      } else if(data->errval[i] < eJK) {
	JK = data->obsval[i]; eJK = data->errval[i]; isJK = 1;
      }
    }
    else if(data->obstype[i] == OBS_TYPE_HKCOLOR) {
      if(!isHK) { 
	HK = data->obsval[i]; eHK = data->errval[i]; isHK = 1;
	Ncolor++;
      } else if(data->errval[i] < eHK) {
	HK = data->obsval[i]; eHK = data->errval[i]; isHK = 1;
      }
    }
  }

  if(Nmag > 0 && Ncolor > 0) {
    for(i=0; i < 4; i++) {
      if(isV) {
	if(!isJ && isVJ) {
	  Jmag = Vmag - VJ;
	  eJmag = sqrt(eVmag*eVmag + eVJ*eVJ);
	  isJ = 1;
	  Nmag++;
	}
	if(!isH && isVH) {
	  Hmag = Vmag - VH;
	  eHmag = sqrt(eHmag*eHmag +eVH*eVH);
	  isH = 1;
	  Nmag++;
	}
	if(!isK && isVK) {
	  Kmag = Vmag - VK;
	  eKmag = sqrt(eKmag*eKmag + eVK*eVK);
	  isK = 1;
	  Nmag++;
	}
      }
      if(isJ) {
	if(!isV && isVJ) {
	  Vmag = Jmag + VJ;
	  eVmag = sqrt(eJmag*eJmag + eVJ*eVJ);
	  isV = 1;
	  Nmag++;
	}
	if(!isH && isJH) {
	  Hmag = Jmag - JH;
	  eHmag = sqrt(eJmag*eJmag + eJH*eJH);
	  isH = 1;
	  Nmag++;
	}
	if(!isK && isJK) {
	  Kmag = Jmag - JK;
	  eKmag = sqrt(eJmag*eJmag + eJK*eJK);
	  isK = 1;
	  Nmag++;
	}
      }
      if(isH) {
	if(!isV && isVH) {
	  Vmag = Hmag + VH;
	  eVmag = sqrt(eHmag*eHmag + eVH*eVH);
	  isV = 1;
	  Nmag++;
	}
	if(!isJ && isJH) {
	  Jmag = Hmag + JH;
	  eJmag = sqrt(eHmag*eHmag + eJH*eJH);
	  isJ = 1;
	  Nmag++;
	}
	if(!isK && isHK) {
	  Kmag = Hmag - HK;
	  eKmag = sqrt(eHmag*eHmag + eHK*eHK);
	  isK = 1;
	  Nmag++;
	}
      }
      if(isK) {
	if(!isV && isVK) {
	  Vmag = Kmag + VK;
	  eVmag = sqrt(eKmag*eKmag + eVK*eVK);
	  isV = 1;
	  Nmag++;
	}
	if(!isJ && isJK) {
	  Jmag = Kmag + JK;
	  eJmag = sqrt(eKmag*eKmag + eJK*eJK);
	  isJ = 1;
	  Nmag++;
	}
	if(!isH && isHK) {
	  Hmag = Kmag + HK;
	  eHmag = sqrt(eHmag*eHmag + eHK*eHK);
	  isH = 1;
	  Nmag++;
	}
      }
    }
  }

  if(!founddist && isAV && Nmag > 1) {
    InitializeEmpiricalModelParams(&mod, 0);
    EvaluateEmpiricalModel(&mdata, &mod, logrhoinit, 0., 0.);
    if(isV) {
      a1_ += (Vmag - mdata.absVmag)/eVmag/eVmag;
      b1_ += 1./eVmag/eVmag;
      c1_ += 1./eVmag/eVmag;
      a2_ += (Vmag - mdata.absVmag)/eVmag/eVmag;
      b2_ += 1./eVmag/eVmag;
      c2_ += 1./eVmag/eVmag;
    }
    if(isJ) {
      a1_ += (Jmag - mdata.absJmag)*REDLAW_AJ/eJmag/eJmag;
      b1_ += (REDLAW_AJ*REDLAW_AJ)/eJmag/eJmag;
      c1_ += (REDLAW_AJ)/eJmag/eJmag;
      a2_ += (Jmag - mdata.absJmag)/eJmag/eJmag;
      b2_ += (REDLAW_AJ)/eJmag/eJmag;
      c2_ += 1./eJmag/eJmag;
    }
    if(isH) {
      a1_ += (Hmag - mdata.absHmag)*REDLAW_AJ/eHmag/eHmag;
      b1_ += (REDLAW_AH*REDLAW_AH)/eHmag/eHmag;
      c1_ += (REDLAW_AH)/eHmag/eHmag;
      a2_ += (Hmag - mdata.absHmag)/eHmag/eHmag;
      b2_ += (REDLAW_AH)/eHmag/eHmag;
      c2_ += 1./eHmag/eHmag;
    }
    if(isK) {
      a1_ += (Kmag - mdata.absKmag)*REDLAW_AK/eKmag/eKmag;
      b1_ += (REDLAW_AK*REDLAW_AK)/eKmag/eKmag;
      c1_ += (REDLAW_AK)/eKmag/eKmag;
      a2_ += (Kmag - mdata.absKmag)/eKmag/eKmag;
      b2_ += (REDLAW_AK)/eKmag/eKmag;
      c2_ += 1./eKmag/eKmag;
    }
    AV = (c1_*a2_ - c2_*a1_)/(c1_*b2_ - b1_*c2_);
    if(AV >= 0.0 && AV <= maxAV) {
      distmod = (a1_ - AV*b1_)/c1_;
      eAV = sqrt(c2_/(c2_*b1_ - c1_*b2_));
      edistmod = sqrt(b1_/(c2_*b1_ - c1_*b2_));
    } else if(AV > maxAV) {
      AV = maxAV;
      a2_ = 0.; c2_ = 0.;
      if(isV) {
	a2_ += (Vmag - mdata.absVmag - AV)/eVmag/eVmag;
	c2_ += 1./eVmag/eVmag;
      }
      if(isJ) {
	a2_ += (Jmag - mdata.absJmag - AV*REDLAW_AJ)/eJmag/eJmag;
	c2_ += 1./eJmag/eJmag;
      }
      if(isH) {
	a2_ += (Hmag - mdata.absHmag - AV*REDLAW_AH)/eHmag/eHmag;
	c2_ += 1./eHmag/eHmag;
      }
      if(isK) {
	a2_ += (Kmag - mdata.absKmag - AV*REDLAW_AK)/eKmag/eKmag;
	c2_ += 1./eKmag/eKmag;
      }
      eAV = 0.1*maxAV;
      distmod = (a2_/c2_);
      edistmod = sqrt(1./c2_);
    } else {
      AV = 0.0;
      eAV = 0.1*maxAV;
      distmod = (a2_/c2_);
      edistmod = sqrt(1./c2_);
    }
    *initdistmod = distmod; *initdistmoderr = edistmod;
    *initAV = AV; *initAVerr = eAV;
    return;
  }
  if(!founddist && isAV && Nmag == 1 && Ncolor > 0) {
    /* First get the reddenning estimate from the color */
    InitializeEmpiricalModelParams(&mod, 0);
    EvaluateEmpiricalModel(&mdata, &mod, logrhoinit, 0., 0.);
    a2_ = 0.; c2_ = 0.;
    if(isVJ) {
      a2_ += (VJ - (mdata.absVmag - mdata.absJmag))*(1. - REDLAW_AJ)/eVJ/eVJ;
      c2_ += (1. - REDLAW_AJ)*(1. - REDLAW_AJ)/eVJ/eVJ;
    }
    if(isVH) {
      a2_ += (VH - (mdata.absVmag - mdata.absHmag))*(1. - REDLAW_AH)/eVH/eVH;
      c2_ += (1. - REDLAW_AH)*(1. - REDLAW_AH)/eVH/eVH;
    }
    if(isVK) {
      a2_ += (VK - (mdata.absVmag - mdata.absKmag))*(1. - REDLAW_AK)/eVK/eVK;
      c2_ += (1. - REDLAW_AK)*(1. - REDLAW_AK)/eVK/eVK;
    }
    if(isJH) {
      a2_ += (JH - (mdata.absJmag - mdata.absHmag))*(REDLAW_AJ - REDLAW_AH)/eJH/eJH;
      c2_ += (REDLAW_AJ - REDLAW_AH)*(REDLAW_AJ - REDLAW_AH)/eJH/eJH;
    }
    if(isJK) {
      a2_ += (JK - (mdata.absJmag - mdata.absKmag))*(REDLAW_AJ - REDLAW_AK)/eJK/eJK;
      c2_ += (REDLAW_AJ - REDLAW_AK)*(REDLAW_AJ - REDLAW_AK)/eJK/eJK;
    }
    if(isHK) {
      a2_ += (HK - (mdata.absHmag - mdata.absKmag))*(REDLAW_AH - REDLAW_AK)/eHK/eHK;
      c2_ += (REDLAW_AH - REDLAW_AK)*(REDLAW_AH - REDLAW_AK)/eHK/eHK;
    }
    AV = (a2_/c2_);
    eAV = sqrt(1./c2_);

    if(AV < 0) { AV = 0.; if(eAV <= 0) {eAV = 0.1*maxAV;} if(eAV > maxAV) {eAV = 0.5*maxAV;}}
    if(AV > maxAV) {AV = maxAV; if(eAV > maxAV) {eAV = 0.5*maxAV;}}

    a2_ = 0.; c2_ = 0.;
    if(isV) {
      a2_ += (Vmag - mdata.absVmag - AV)/eVmag/eVmag;
      c2_ += 1./eVmag/eVmag;
    }
    if(isJ) {
      a2_ += (Jmag - mdata.absJmag - AV*REDLAW_AJ)/eJmag/eJmag;
      c2_ += 1./eJmag/eJmag;
    }
    if(isH) {
      a2_ += (Hmag - mdata.absHmag - AV*REDLAW_AH)/eHmag/eHmag;
      c2_ += 1./eHmag/eHmag;
    }
    if(isK) {
      a2_ += (Kmag - mdata.absKmag - AV*REDLAW_AK)/eKmag/eKmag;
      c2_ += 1./eKmag/eKmag;
    }
    distmod = (a2_/c2_);
    edistmod = sqrt(1./c2_);
    
    *initdistmod = distmod; *initdistmoderr = edistmod;
    *initAV = AV; *initAVerr = eAV;
    return;
  }    

  if(isAV && Nmag == 0 && Ncolor > 0) {
    /* Not enough information for the distance. But we can get the reddenning estimate from the color */
    InitializeEmpiricalModelParams(&mod, 0);
    EvaluateEmpiricalModel(&mdata, &mod, logrhoinit, 0., 0.);
    a2_ = 0.; c2_ = 0.;
    if(isVJ) {
      a2_ += (VJ - (mdata.absVmag - mdata.absJmag))*(1. - REDLAW_AJ)/eVJ/eVJ;
      c2_ += (1. - REDLAW_AJ)*(1. - REDLAW_AJ)/eVJ/eVJ;
    }
    if(isVH) {
      a2_ += (VH - (mdata.absVmag - mdata.absHmag))*(1. - REDLAW_AH)/eVH/eVH;
      c2_ += (1. - REDLAW_AH)*(1. - REDLAW_AH)/eVH/eVH;
    }
    if(isVK) {
      a2_ += (VK - (mdata.absVmag - mdata.absKmag))*(1. - REDLAW_AK)/eVK/eVK;
      c2_ += (1. - REDLAW_AK)*(1. - REDLAW_AK)/eVK/eVK;
    }
    if(isJH) {
      a2_ += (JH - (mdata.absJmag - mdata.absHmag))*(REDLAW_AJ - REDLAW_AH)/eJH/eJH;
      c2_ += (REDLAW_AJ - REDLAW_AH)*(REDLAW_AJ - REDLAW_AH)/eJH/eJH;
    }
    if(isJK) {
      a2_ += (JK - (mdata.absJmag - mdata.absKmag))*(REDLAW_AJ - REDLAW_AK)/eJK/eJK;
      c2_ += (REDLAW_AJ - REDLAW_AK)*(REDLAW_AJ - REDLAW_AK)/eJK/eJK;
    }
    if(isHK) {
      a2_ += (HK - (mdata.absHmag - mdata.absKmag))*(REDLAW_AH - REDLAW_AK)/eHK/eHK;
      c2_ += (REDLAW_AH - REDLAW_AK)*(REDLAW_AH - REDLAW_AK)/eHK/eHK;
    }
    AV = (a2_/c2_);
    eAV = sqrt(1./c2_);

    if(AV < 0) { AV = 0.; if(eAV <= 0) {eAV = 0.1*maxAV;} if(eAV > maxAV) {eAV = 0.5*maxAV;}}
    if(AV > maxAV) {AV = maxAV; if(eAV > maxAV) {eAV = 0.5*maxAV;}}

    distmod = 0.;
    edistmod = 5.;
    
    *initdistmod = distmod; *initdistmoderr = edistmod;
    *initAV = AV; *initAVerr = eAV;
    return;
  }    

  if(!founddist && !isAV && Nmag > 0) {
    /* No reddening, just compute the distance modulus from the magnitudes */
    InitializeEmpiricalModelParams(&mod, 0);
    EvaluateEmpiricalModel(&mdata, &mod, logrhoinit, 0., 0.);
    a2_ = 0.; c2_ = 0.;
    if(isV) {
      a2_ += (Vmag - mdata.absVmag)/eVmag/eVmag;
      c2_ += 1./eVmag/eVmag;
    }
    if(isJ) {
      a2_ += (Jmag - mdata.absJmag)/eJmag/eJmag;
      c2_ += 1./eJmag/eJmag;
    }
    if(isH) {
      a2_ += (Hmag - mdata.absHmag)/eHmag/eHmag;
      c2_ += 1./eHmag/eHmag;
    }
    if(isK) {
      a2_ += (Kmag - mdata.absKmag)/eKmag/eKmag;
      c2_ += 1./eKmag/eKmag;
    }
    AV = 0.;
    eAV = 0.;
    distmod = (a2_/c2_);
    edistmod = sqrt(1./c2_);
    *initdistmod = distmod; *initdistmoderr = edistmod;
    *initAV = AV; *initAVerr = eAV;
    return;
  }

  if(founddist && isAV && Nmag > 0 && Ncolor == 0) {
    /* distance modulus constrained as an observable, compute AV from the
       magnitudes */
    InitializeEmpiricalModelParams(&mod, 0);
    EvaluateEmpiricalModel(&mdata, &mod, logrhoinit, 0., 0.);
    if(isV) {
      a2_ += (Vmag - (mdata.absVmag + distmod))/(eVmag*eVmag + edistmod*edistmod);
      c2_ += 1./(eVmag*eVmag + edistmod*edistmod);
    }
    if(isJ) {
      a2_ += (Jmag - (mdata.absJmag + distmod))*REDLAW_AJ/(eJmag*eJmag + edistmod*edistmod);
      c2_ += REDLAW_AJ*REDLAW_AJ/(eJmag*eJmag + edistmod*edistmod);
    }
    if(isH) {
      a2_ += (Hmag - (mdata.absHmag + distmod))*REDLAW_AH/(eHmag*eHmag + edistmod*edistmod);
      c2_ += REDLAW_AH*REDLAW_AH/(eHmag*eHmag + edistmod*edistmod);
    }
    if(isK) {
      a2_ += (Kmag - (mdata.absKmag + distmod))*REDLAW_AK/(eKmag*eKmag + edistmod*edistmod);
      c2_ += REDLAW_AK*REDLAW_AK/(eKmag*eKmag + edistmod*edistmod);
    }
    AV = a2_/c2_;
    eAV = sqrt(1./c2_);
    if(AV <= 0.0) AV = 0.;
    if(AV > maxAV) {AV = maxAV; eAV = 0.5*maxAV;}
    *initdistmod = distmod; *initdistmoderr = edistmod;
    *initAV = AV; *initAVerr = eAV;
    return;
  }
  if(!founddist) {
    *initdistmod = 0.;
    *initdistmoderr = 5.;
  }
  *initAV = 0.; *initAVerr = 0.1;
  return;
}


/* This function determines the initial estimates for the density */
void GuessInitialDensity(double *initlogrho, double *initlogrhoerr, _ObsData *data) {
  int i;

  int distindex, magindex, magtype, magindex2, angdiamindex;
  double distmod, edistmod;

  /* Determine initial estimates for logrho based on any observable parameters */
  for(i=0; i < data->Nobs; i++) {
    if(data->obstype[i] == OBS_TYPE_DENSITY) {
      *initlogrho = log(data->obsval[i])/M_LN10;
      if(data->obsval[i] > 0.0) {
	*initlogrhoerr = data->errval[i]/data->obsval[i]/M_LN10;
      } else {
	*initlogrhoerr = 0.1;
      }
      return;
    }
    else if(data->obstype[i] == OBS_TYPE_LOGDENSITY) {
      *initlogrho = data->obsval[i];
      *initlogrhoerr = data->errval[i];
      return;
    }
    else if(data->obstype[i] == OBS_TYPE_DENSITYSOLAR) {
      *initlogrho = log(data->obsval[i]*SOLARDENSITY)/M_LN10;
      if(data->obsval[i] > 0.0) {
	*initlogrhoerr = data->errval[i]/data->obsval[i]/M_LN10;
      } else {
	*initlogrhoerr = 0.1;
      }
      return;
    }
    else if(data->obstype[i] == OBS_TYPE_LOGDENSITYSOLAR) {
      *initlogrho = data->obsval[i] + log(SOLARDENSITY)/M_LN10;
      *initlogrhoerr = data->errval[i];
      return;
    }
  }
  /* The density isn't an initial parameter; check for radius. */
  for(i=0; i < data->Nobs; i++) {
    if(data->obstype[i] == OBS_TYPE_RADIUS) {
      GetDensityFromValue(initlogrho, initlogrhoerr,data->obstype[i], data->obsval[i], data->errval[i]);
      return;
    }
  }
  /* Check for mass */
  for(i=0; i < data->Nobs; i++) {
    if(data->obstype[i] == OBS_TYPE_MASS) {
      GetDensityFromValue(initlogrho, initlogrhoerr, OBS_TYPE_MASS, data->obsval[i], data->errval[i]);
      return;
    }
  }
  /* Check for absolute magnitudes J, H or K */
  for(i=0; i < data->Nobs; i++) {
    if(data->obstype[i] == OBS_TYPE_ABSJMAG || 
       data->obstype[i] == OBS_TYPE_ABSHMAG ||
       data->obstype[i] == OBS_TYPE_ABSKMAG) {
      GetDensityFromValue(initlogrho, initlogrhoerr, data->obstype[i], data->obsval[i], data->errval[i]);
      return;
    }
  }
  /* Check for absolute V magnitude */
  for(i=0; i < data->Nobs; i++) {
    if(data->obstype[i] == OBS_TYPE_ABSVMAG) {
      GetDensityFromValue(initlogrho, initlogrhoerr, data->obstype[i], data->obsval[i], data->errval[i]);
      return;
    }
  }

  /* Check for distance plus a magnitude */
  distindex = -1; magindex = -1; angdiamindex = -1;
  for(i=0; i < data->Nobs; i++) {
    if(data->obstype[i] == OBS_TYPE_KMAG ||
       data->obstype[i] == OBS_TYPE_JMAG ||
       data->obstype[i] == OBS_TYPE_HMAG) {
      magindex = i;
      if(data->obstype[i] == OBS_TYPE_KMAG)
	magtype = OBS_TYPE_ABSKMAG;
      else if(data->obstype[i] == OBS_TYPE_JMAG) 
	magtype = OBS_TYPE_ABSJMAG;
      else if(data->obstype[i] == OBS_TYPE_HMAG)
	magtype = OBS_TYPE_ABSHMAG;
    }
    else if(data->obstype[i] == OBS_TYPE_VMAG &&
	    magindex < 0) {
      magindex = i;
      magtype = OBS_TYPE_ABSVMAG;
    }
    else if(data->obstype[i] == OBS_TYPE_DIST) {
      distindex = i;
      distmod = 5*log(data->obsval[i])/M_LN10 - 5.0;
      edistmod = 5*data->errval[i]/data->obsval[i]/M_LN10;
    }
    else if(data->obstype[i] == OBS_TYPE_DISTMOD) {
      distindex = i;
      distmod = data->obsval[i];
      edistmod = data->errval[i];
    }
    else if(data->obstype[i] == OBS_TYPE_PARALLAX) {
      distindex = i;
      distmod = 5*log((1000.0/data->obsval[i]))/M_LN10 - 5.0;
      edistmod = 5*data->errval[i]/data->obsval[i]/M_LN10;
    }
    else if(data->obstype[i] == OBS_TYPE_ANGDIAM) {
      angdiamindex = i;
    }
  }
  if(magindex >= 0 && distindex >= 0) {
    GetDensityFromValue(initlogrho, initlogrhoerr, magtype, (data->obsval[magindex]-distmod), (sqrt(data->errval[magindex]*data->errval[magindex] + edistmod*edistmod)));
    return;
  }
  else if(distindex >= 0 && angdiamindex >= 0) {
    GetDensityFromValue(initlogrho, initlogrhoerr, OBS_TYPE_RADIUS, (data->obsval[angdiamindex]*exp(M_LN10*((distmod + 5.0)/5.0))/9.298438), (sqrt(data->errval[angdiamindex]*data->errval[angdiamindex]*exp(2*M_LN10*((distmod + 5.0)/5.0))/(9.298438*9.298438) + edistmod*edistmod*data->obsval[angdiamindex]*data->obsval[angdiamindex]*M_LN10*M_LN10*exp(2*M_LN10*((distmod+5.0)/5.0))/5.0/5.0/9.298438/9.298438)));
    return;
  }
  
  /* Check for Teff */
  for(i=0; i < data->Nobs; i++) {
    if(data->obstype[i] == OBS_TYPE_TEFF) {
      GetDensityFromValue(initlogrho, initlogrhoerr, data->obstype[i], data->obsval[i], data->errval[i]);
      return;
    }
  }

  /* Check for V-K color */ 
  for(i=0; i < data->Nobs; i++) {
    if(data->obstype[i] == OBS_TYPE_VKCOLOR) {
      GetDensityFromValue(initlogrho, initlogrhoerr, data->obstype[i], data->obsval[i], data->errval[i]);
      return;
    }
  }

  /* Check for V and K apparent magnitudes */
  magindex = -1; magindex2 = -1;
  for(i=0; i < data->Nobs; i++) {
    if(data->obstype[i] == OBS_TYPE_VMAG) {
      magindex = i;
    }
    else if(data->obstype[i] == OBS_TYPE_KMAG) {
      magindex2 = i;
    }
  }
  if(magindex >= 0 && magindex2 >= 0) {
    GetDensityFromValue(initlogrho, initlogrhoerr, OBS_TYPE_VKCOLOR, (data->obsval[magindex]-data->obsval[magindex2]), (sqrt(data->errval[magindex]*data->errval[magindex] + data->errval[magindex2]*data->errval[magindex2])));
    return;
  }

  /* Check for V-H color */ 
  for(i=0; i < data->Nobs; i++) {
    if(data->obstype[i] == OBS_TYPE_VHCOLOR) {
      GetDensityFromValue(initlogrho, initlogrhoerr, data->obstype[i], data->obsval[i], data->errval[i]);
      return;
    }
  }

  /* Check for V and H apparent magnitudes */
  magindex = -1; magindex2 = -1;
  for(i=0; i < data->Nobs; i++) {
    if(data->obstype[i] == OBS_TYPE_VMAG) {
      magindex = i;
    }
    else if(data->obstype[i] == OBS_TYPE_HMAG) {
      magindex2 = i;
    }
  }
  if(magindex >= 0 && magindex2 >= 0) {
    GetDensityFromValue(initlogrho, initlogrhoerr, OBS_TYPE_VHCOLOR, (data->obsval[magindex]-data->obsval[magindex2]), (sqrt(data->errval[magindex]*data->errval[magindex] + data->errval[magindex2]*data->errval[magindex2])));
    return;
  }

  /* Check for V-J color */ 
  for(i=0; i < data->Nobs; i++) {
    if(data->obstype[i] == OBS_TYPE_VJCOLOR) {
      GetDensityFromValue(initlogrho, initlogrhoerr, data->obstype[i], data->obsval[i], data->errval[i]);
      return;
    }
  }

  /* Check for V and J apparent magnitudes */
  magindex = -1; magindex2 = -1;
  for(i=0; i < data->Nobs; i++) {
    if(data->obstype[i] == OBS_TYPE_VMAG) {
      magindex = i;
    }
    else if(data->obstype[i] == OBS_TYPE_JMAG) {
      magindex2 = i;
    }
  }
  if(magindex >= 0 && magindex2 >= 0) {
    GetDensityFromValue(initlogrho, initlogrhoerr, OBS_TYPE_VJCOLOR, (data->obsval[magindex]-data->obsval[magindex2]), (sqrt(data->errval[magindex]*data->errval[magindex] + data->errval[magindex2]*data->errval[magindex2])));
    return;
  }

  /* Check for LOGG or LOGGSOLAR */
  for(i=0; i < data->Nobs; i++) {
    if(data->obstype[i] == OBS_TYPE_LOGG ||
       data->obstype[i] == OBS_TYPE_LOGGSOLAR) {
      GetDensityFromValue(initlogrho, initlogrhoerr, data->obstype[i],
			  data->obsval[i], data->errval[i]);
      return;
    }
  }

  /* Check for J-K color */ 
  for(i=0; i < data->Nobs; i++) {
    if(data->obstype[i] == OBS_TYPE_JKCOLOR) {
      GetDensityFromValue(initlogrho, initlogrhoerr, data->obstype[i], data->obsval[i], data->errval[i]);
      return;
    }
  }

  /* Check for J and K apparent magnitudes */
  magindex = -1; magindex2 = -1;
  for(i=0; i < data->Nobs; i++) {
    if(data->obstype[i] == OBS_TYPE_JMAG) {
      magindex = i;
    }
    else if(data->obstype[i] == OBS_TYPE_KMAG) {
      magindex2 = i;
    }
  }
  if(magindex >= 0 && magindex2 >= 0) {
    GetDensityFromValue(initlogrho, initlogrhoerr, OBS_TYPE_JKCOLOR, (data->obsval[magindex]-data->obsval[magindex2]), (sqrt(data->errval[magindex]*data->errval[magindex] + data->errval[magindex2]*data->errval[magindex2])));
    return;
  }

  /* Check for J-H color */ 
  for(i=0; i < data->Nobs; i++) {
    if(data->obstype[i] == OBS_TYPE_JHCOLOR) {
      GetDensityFromValue(initlogrho, initlogrhoerr, data->obstype[i], data->obsval[i], data->errval[i]);
      return;
    }
  }

  /* Check for J and H apparent magnitudes */
  magindex = -1; magindex2 = -1;
  for(i=0; i < data->Nobs; i++) {
    if(data->obstype[i] == OBS_TYPE_JMAG) {
      magindex = i;
    }
    else if(data->obstype[i] == OBS_TYPE_HMAG) {
      magindex2 = i;
    }
  }
  if(magindex >= 0 && magindex2 >= 0) {
    GetDensityFromValue(initlogrho, initlogrhoerr, OBS_TYPE_JHCOLOR, (data->obsval[magindex]-data->obsval[magindex2]), (sqrt(data->errval[magindex]*data->errval[magindex] + data->errval[magindex2]*data->errval[magindex2])));
    return;
  }

  /* Check for H-K color */ 
  for(i=0; i < data->Nobs; i++) {
    if(data->obstype[i] == OBS_TYPE_HKCOLOR) {
      GetDensityFromValue(initlogrho, initlogrhoerr, data->obstype[i], data->obsval[i], data->errval[i]);
      return;
    }
  }

  /* Check for H and K apparent magnitudes */
  magindex = -1; magindex2 = -1;
  for(i=0; i < data->Nobs; i++) {
    if(data->obstype[i] == OBS_TYPE_HMAG) {
      magindex = i;
    }
    else if(data->obstype[i] == OBS_TYPE_KMAG) {
      magindex2 = i;
    }
  }
  if(magindex >= 0 && magindex2 >= 0) {
    GetDensityFromValue(initlogrho, initlogrhoerr, OBS_TYPE_HKCOLOR, (data->obsval[magindex]-data->obsval[magindex2]), (sqrt(data->errval[magindex]*data->errval[magindex] + data->errval[magindex2]*data->errval[magindex2])));
    return;
  }

  /* No good density constraint from the observables given */
  *initlogrho = 1.0;
  *initlogrhoerr = 0.75;

  return;
}

void CopyModelData(_ModelData *outmodel,_ModelData *inmodel) {
  outmodel->logrho = inmodel->logrho;
  outmodel->dlR = inmodel->dlR;
  outmodel->dMbol = inmodel->dMbol;
  outmodel->lM = inmodel->lM;
  outmodel->lR = inmodel->lR;
  outmodel->lTeff = inmodel->lTeff;
  outmodel->Mbol = inmodel->Mbol;
  outmodel->BCV = inmodel->BCV;
  outmodel->BCJ = inmodel->BCJ;
  outmodel->BCH = inmodel->BCH;
  outmodel->BCK = inmodel->BCK;
  outmodel->absVmag = inmodel->absVmag;
  outmodel->absJmag = inmodel->absJmag;
  outmodel->absHmag = inmodel->absHmag;
  outmodel->absKmag = inmodel->absKmag;
  outmodel->AV = inmodel->AV;
  outmodel->distmod = inmodel->distmod;
}


/* Run the differential evolution markov chain monte carlo fit for the stellar
   parameters.  The input is the structure of observed data.
   Output will be the MCMC chain in outchains, the likelihoods associated
   with each link in the chain in outloglike, and the number of links in
   Nlinks. */
void RunDEMCMC(_ObsData *data, _ModelData **outchains, double **outloglike, int *Nlinks, int useAV, double maxAV, int nlinkstest, int maxlinks, double maxmem, int *outisdistmod, int *outisAV) {
  int nchain = 0;
  int nvar = 0;
  int i, isAV, isdistmod, j;
  int NGR_pass = 0;
  int ngr_cnt, ngr_cnt_last;
  double gamma, eps, initlogrho, initlogrhoerr, alpha, u;
  double W;
  double B;
  double R;
  double theta_avg;
  int test_pass;

  double **bvector; /* vector of unknown variables; 1 per chain */
  double **thetavector; /* vector for use in Gelman & Rubin statistic */
  double **s2vector; /* vector of s2 values, to use in calculating 
                        Gelman & Rubin statistic */
  double *chi2vector; /* -2lnL values for each chain */
  double *dvector;
  double *gvector;
  double *errvector;
  _ModelData *mdata, testmdata; /* Stellar models for each chain */

  _EmpiricalModelParams testmod;

  int size_out_chains = 0, k, ntot, nindx, cur_chain_length, ichain, ccnt, is_not_accepted, jchain, kchain;
  int maxNalloc;

  int iiter;

  double initdistmod, initdistmoderr, initAV, initAVerr, nchi, vx;
  double tmp1, tmp2;

  double **simplex = NULL;
  double *simplex_chi2 = NULL;
  int Nparam_simplex = 0;
  int Ntovary_simplex =0;
  int *ia_simplex = NULL;
  int nfunkeval;

  double chi2best, *red, redmed;
  int ibest;

  int amoeba_val;

  nvar = 3 + 22;
  isAV = 0; isdistmod = 0;

  /* Check if we need AV or distmod as free parameters */
  for(i=0; i < data->Nobs; i++) {
    if(data->obstype[i] == OBS_TYPE_VMAG ||
       data->obstype[i] == OBS_TYPE_JMAG ||
       data->obstype[i] == OBS_TYPE_HMAG ||
       data->obstype[i] == OBS_TYPE_KMAG ||
       data->obstype[i] == OBS_TYPE_DISTMOD ||
       data->obstype[i] == OBS_TYPE_DIST ||
       data->obstype[i] == OBS_TYPE_PARALLAX ||
       data->obstype[i] == OBS_TYPE_ANGDIAM) {
      isdistmod = 1;
      if(useAV) isAV = 1;
    }
    if(data->obstype[i] == OBS_TYPE_AV || data->obstype[i] == OBS_TYPE_EBV) {
      isAV = 1;
    }
    if(data->obstype[i] == OBS_TYPE_VJCOLOR ||
       data->obstype[i] == OBS_TYPE_VHCOLOR ||
       data->obstype[i] == OBS_TYPE_VKCOLOR ||
       data->obstype[i] == OBS_TYPE_JHCOLOR ||
       data->obstype[i] == OBS_TYPE_JKCOLOR ||
       data->obstype[i] == OBS_TYPE_HKCOLOR) {
      if(useAV) isAV = 1;
    }
  }
  if(isAV) nvar++;
  if(isdistmod) nvar++;

  *outisdistmod = isdistmod;
  *outisAV = isAV;

  nchain = 2*nvar;
  if(nchain < 3) nchain = 3;

  gamma = 1.7/sqrt(2.0*nvar);
  eps = 0.0001;

  size_out_chains = 2*ceil(nlinkstest/2.);

  maxNalloc = floor((maxmem*1.0e9)/(sizeof(_ModelData)+sizeof(double)));

  if(size_out_chains > maxNalloc) size_out_chains = maxNalloc;

  if(((*outchains) = (_ModelData *) malloc(size_out_chains*sizeof(_ModelData))) == NULL ||
     ((*outloglike) = (double *) malloc(size_out_chains*sizeof(double))) == NULL)
    {
      fprintf(stderr,"Memory Allocation Error\n");
      exit(1);
    }

  /* bvector is the vector of unkown variables; 1 vector per chain is
     required. */
  if((bvector = (double **) malloc(nchain * sizeof(double *))) == NULL ||
     (s2vector = (double **) malloc(nchain * sizeof(double *))) == NULL ||
     (chi2vector = (double *) malloc(nchain * sizeof(double))) == NULL ||
     (dvector = (double *) malloc(nvar * sizeof(double))) == NULL ||
     (gvector = (double *) malloc(nvar * sizeof(double))) == NULL ||
     (mdata = (_ModelData *) malloc(nchain * sizeof(_ModelData))) == NULL ||
     (thetavector = (double **) malloc(nchain * sizeof(double *))) == NULL ||
     (errvector = (double *) malloc(nvar * sizeof(double))) == NULL ||
     (red = (double *) malloc(nchain * sizeof(double))) == NULL) {
    fprintf(stderr,"Memory Allocation Error\n");
    exit(1);
  }
  for(i=0; i < nchain; i++) {
    if((bvector[i] = (double *) malloc(nvar * sizeof(double))) == NULL ||
       (thetavector[i] = (double *) malloc(nvar * sizeof(double))) == NULL ||
       (s2vector[i] = (double *) malloc(nvar * sizeof(double))) == NULL) {
      fprintf(stderr,"Memory Allocation Error\n");
    }
  }

  InitializeEmpiricalModelParams(&testmod, 0);
  GuessInitialDensity(&initlogrho,&initlogrhoerr,data);

  /* First run a DHSX fit */

  incrementparameters_foramoeba(&Nparam_simplex, &Ntovary_simplex, &simplex, &ia_simplex, 1, initlogrho, initlogrhoerr);
  incrementparameters_foramoeba(&Nparam_simplex, &Ntovary_simplex, &simplex, &ia_simplex, 1, 0., testmod.v[INDEX_Sr]);
  incrementparameters_foramoeba(&Nparam_simplex, &Ntovary_simplex, &simplex, &ia_simplex, 1, 0., testmod.v[INDEX_Sm]);
  for(j=3; j < NEMPIRICALMODELPARAMS+3; j++) {
    incrementparameters_foramoeba(&Nparam_simplex, &Ntovary_simplex, &simplex, &ia_simplex, 1, testmod.v[j-3], ModelStddev[j-3]);
  }
  if(isdistmod || isAV) {
    GuessInitialDistAV(&initdistmod,&initdistmoderr,&initAV,&initAVerr,isdistmod,isAV,data,maxAV,initlogrho);
  }
  if(isdistmod) {
    incrementparameters_foramoeba(&Nparam_simplex, &Ntovary_simplex, &simplex, &ia_simplex, 1, initdistmod, initdistmoderr);
  }
  if(isAV) {
    incrementparameters_foramoeba(&Nparam_simplex, &Ntovary_simplex, &simplex, &ia_simplex, 1, initAV, initAVerr);
  }
  
  amoeba_initializesimplexchi2(Nparam_simplex, Ntovary_simplex, simplex, &simplex_chi2, data, &(mdata[0]), isAV, isdistmod, maxAV);

  amoeba_val = amoeba(simplex, simplex_chi2, ia_simplex, Nparam_simplex, 1e-10, &nfunkeval, 100000, data, &(mdata[0]), isAV, isdistmod, maxAV);

  if(!amoeba_val) {
    chi2best = simplex_chi2[0];
    ibest = 0;
    for(i=1; i < Ntovary_simplex + 1; i++) {
      if(!(simplex_chi2[i] < BADCHI2) && (simplex_chi2[i] < chi2best || (chi2best < BADCHI2))) {
	chi2best = simplex_chi2[i];
	ibest = i;
      }
    }

    /* Initialize the bvectors */
    j = 0;
    bvector[0][0] = simplex[ibest][0]; /*initlogrho;*/
    bvector[0][1] = simplex[ibest][1]; /*0.;*/  /* dlR */
    bvector[0][2] = simplex[ibest][2]; /*0.;*/  /* dlM */
    for(j=3; j < NEMPIRICALMODELPARAMS+3; j++) {
      bvector[0][j] = simplex[ibest][j]; /*testmod.v[j-3];*/
    }
    /* distmod and AV if they are free parameters */
    /*  if(isdistmod || isAV) {
	GuessInitialDistAV(&initdistmod,&initdistmoderr,&initAV,&initAVerr,isdistmod,isAV,data,maxAV,initlogrho);
	}*/
    if(isdistmod) {
      bvector[0][j] = simplex[ibest][j]; /*initdistmod;*/
      j++;
    }
    if(isAV) {
      bvector[0][j] = simplex[ibest][j]; /*initAV;*/
      j++;
    }
    
    chi2vector[0] = GetLnLikelihood(data, bvector[0], &(mdata[0]), isAV, isdistmod, maxAV);
    
    if(chi2vector[0] <= BADCHI2) {
      tmp1 = MIN_LRHO_LIMIT;
      tmp2 = MAX_LRHO_LIMIT;
      fprintf(stderr,"#WARNING: initial parameters are outside of the constraints (either supplied limits on one or more of the observations, or else log10(rho) < %.1f or log10(rho) > %.1f.\n", tmp1, tmp2);
    }

    errvector[0] = GetAmoeba_Uncertainties(simplex[ibest], 0, Nparam_simplex, simplex_chi2[ibest], initlogrhoerr, data, &(mdata[0]), isAV, isdistmod, maxAV);
    errvector[1] = GetAmoeba_Uncertainties(simplex[ibest], 1, Nparam_simplex, simplex_chi2[ibest], testmod.v[INDEX_Sr], data, &(mdata[0]), isAV, isdistmod, maxAV);
    errvector[2] = GetAmoeba_Uncertainties(simplex[ibest], 2, Nparam_simplex, simplex_chi2[ibest], testmod.v[INDEX_Sm], data, &(mdata[0]), isAV, isdistmod, maxAV);
    for(j=3; j < NEMPIRICALMODELPARAMS+3; j++) {
      errvector[j] = ModelStddev[j-3];
    }
    if(isdistmod) {errvector[j] = GetAmoeba_Uncertainties(simplex[ibest], j, Nparam_simplex, simplex_chi2[ibest], initdistmoderr, data, &(mdata[0]), isAV, isdistmod, maxAV); j++;}
    if(isAV) {errvector[j] = GetAmoeba_Uncertainties(simplex[ibest], j, Nparam_simplex, simplex_chi2[ibest], initAVerr, data, &(mdata[0]), isAV, isdistmod, maxAV); j++;}
  }
  else {
    do {
      bvector[0][0] = initlogrho + gaussian()*initlogrhoerr;
      bvector[0][1] = gaussian()*testmod.v[INDEX_Sr];  /* dlR */
      bvector[0][2] = gaussian()*testmod.v[INDEX_Sm];  /* dlM */
      for(j=3; j < NEMPIRICALMODELPARAMS+3; j++) {
	bvector[0][j] = testmod.v[j-3];
      }
      if(isdistmod) {
	bvector[0][j] = initdistmod + gaussian()*initdistmoderr;
	j++;
      }
      if(isAV) {
	bvector[0][j] = initAV + gaussian()*initAVerr;
	j++;
      }
      chi2vector[0] = GetLnLikelihood(data, bvector[0], &(mdata[0]), isAV, isdistmod, maxAV);
    } while (chi2vector[0] <= BADCHI2);
    errvector[0] = initlogrhoerr;
    errvector[1] = testmod.v[INDEX_Sr];
    errvector[1] = testmod.v[INDEX_Sm];
    for(j=3; j < NEMPIRICALMODELPARAMS+3; j++) {
      errvector[j] = ModelStddev[j-3];
    }
    if(isdistmod) {errvector[j] = initdistmoderr; j++;}
    if(isAV) {errvector[j] = initAVerr; j++;}
  }

  for(k=1; k < nchain; k++) {
    InitializeEmpiricalModelParams(&testmod, 1);
    red[k] = 1.0;
    do {
      bvector[k][0] = bvector[0][0] + gaussian()*errvector[0]/red[k];
      bvector[k][1] = bvector[0][1] + gaussian()*errvector[1]/red[k];  /* dlR */
      bvector[k][2] = bvector[0][2] + gaussian()*errvector[2]/red[k];  /* dlM */
      for(j=3; j < NEMPIRICALMODELPARAMS+3; j++) {
	bvector[k][j] = bvector[0][j] + gaussian()*(testmod.v[j-3] - bvector[0][j])/red[k]; /*testmod.v[j-3];*/
      }
      if(isdistmod) {
	bvector[k][j] = bvector[0][j] + gaussian()*errvector[j]/red[k];
	j++;
      }
      if(isAV) {
	bvector[k][j] = bvector[0][j] + gaussian()*errvector[j]/red[k];
	j++;
      }
      chi2vector[k] = GetLnLikelihood(data, bvector[k], &(mdata[k]), isAV, isdistmod, maxAV);
      red[k] = red[k] + 1.0;
    } while(chi2vector[k] <= BADCHI2 || (fabs(chi2vector[k] - chi2vector[0]) > 30.0));
  }
  redmed = ComputeMedian(&(red[1]), (nchain-1));
  for(j=0; j < nvar; j++) {
    errvector[j] = errvector[j] / redmed;
  }


  
  /* Now to run the DEMCMC algorithm */

  ntot = 0;
  nindx = 0; 
  cur_chain_length = 0;
  ichain = -1;
  /* Initialize Gelman & Rubin statistics vectors */
  for(i=0; i < nchain; i++) {
    for(j=0; j < nvar; j++) {
      s2vector[i][j] = 0.0;
      thetavector[i][j] = 0.0;
    }
  }
  ngr_cnt = 0;
  while(NGR_pass < NGR_PASS_STOP_LIMIT) {
    if(NGR_pass > 0) {
      ngr_cnt_last = ngr_cnt;
    }
    for(ccnt=0, iiter=0 ; (NGR_pass <= 0 ? ccnt < nlinkstest : iiter < ngr_cnt_last*nchain) ; ntot++, nindx++, iiter++) {
      if(maxlinks > 0 && ntot > maxlinks) break;
      if(nindx >= maxNalloc) nindx = 0;
      ichain++;
      ichain = ichain % nchain;
      is_not_accepted = 0;
    
      /* Choose two random chains to select from */
      jchain = rand() % (nchain - 1);
      kchain = rand() % (nchain - 2);
      if(jchain >= ichain) jchain++;
      if(kchain >= (jchain < ichain ? jchain : ichain)) kchain++;
      if(kchain >= (jchain < ichain ? ichain : jchain)) kchain++;
      
      for(j=0; j < nvar; j++) {
	gvector[j] = gaussian();
      }
      for(i=0; i < nvar; i++) {
	dvector[i] = errvector[i]*gvector[i];
      }
      for(i=0; i < nvar; i++) {
	dvector[i] = dvector[i]*eps;
	dvector[i] = dvector[i] + gamma*(bvector[jchain][i]-bvector[kchain][i]) + bvector[ichain][i];
      }
      
      if(dvector[1+NEMPIRICALMODELPARAMS] < TINYERRTERM) {
	dvector[1+NEMPIRICALMODELPARAMS] = TINYERRTERM;
	is_not_accepted = 1;
      }
      if(dvector[2+NEMPIRICALMODELPARAMS] < TINYERRTERM) {
	dvector[2+NEMPIRICALMODELPARAMS] = TINYERRTERM;
	is_not_accepted = 1;
      }
      if(isAV) {
	if(dvector[nvar-1] < 0) {dvector[nvar-1] = 0.; is_not_accepted = 1;}
	if(dvector[nvar-1] > maxAV) {dvector[nvar-1] = maxAV; is_not_accepted = 1;}
      }
      if(!is_not_accepted) {
	/* Get the likelihood for the trial model */
	nchi = GetLnLikelihood(data, dvector, &testmdata, isAV, isdistmod, maxAV);
	if(nchi > BADCHI2) {
	  vx = exp(-0.5*(nchi-chi2vector[ichain]));
	  alpha=(vx<1.0?vx:1.0);
	  u = drand48();
	}
	if(nchi > BADCHI2 ? (alpha > 0.0 && u <= alpha) : 0) {
	  for(i=0; i < nvar; i++) {
	    bvector[ichain][i] = dvector[i];
	  }
	  chi2vector[ichain]=nchi;
	  CopyModelData(&(mdata[ichain]),&testmdata);
	  
	  ccnt++;
	} else is_not_accepted = 1;
      }
      if(nindx >= size_out_chains) {
	size_out_chains *= 2;
	if(size_out_chains > maxNalloc) size_out_chains = maxNalloc;
	if((*outchains = (_ModelData *) realloc((*outchains),size_out_chains * sizeof(_ModelData))) == NULL ||
	   (*outloglike = (double *) realloc((*outloglike),size_out_chains * sizeof(double))) == NULL) {
	  fprintf(stderr,"Memory Allocation Error\n");
	  exit(1);
	}
      }
      CopyModelData(&((*outchains)[nindx]),&(mdata[ichain]));
      (*outloglike)[nindx] = chi2vector[ichain];
      for(j=0; j < nvar; j++) {
	s2vector[ichain][j] += bvector[ichain][j]*bvector[ichain][j];
	thetavector[ichain][j] += bvector[ichain][j];
      }
      if(ichain == (nchain - 1))
	ngr_cnt++;
      if(cur_chain_length < maxNalloc) cur_chain_length++;
    }
    if(maxlinks > 0 && ntot > maxlinks) break;
    /* Calculate the Gelman & Rubin statistic */
    test_pass = 0;
    if(NGR_pass == 0) {
      for(j=0; j < nvar; j++) {
	W = 0.0;
	B = 0.0;
	theta_avg = 0.0;
	for(i=0; i < nchain; i++) {
	  s2vector[i][j] = s2vector[i][j]/ngr_cnt - (thetavector[i][j]/ngr_cnt)*(thetavector[i][j]/ngr_cnt);
	  W += s2vector[i][j]/nchain;
	  theta_avg += (thetavector[i][j]/ngr_cnt);
	}
	theta_avg = theta_avg/nchain;
	for(i=0; i < nchain; i++) {
	  B += (thetavector[i][j]/ngr_cnt - theta_avg)*(thetavector[i][j]/ngr_cnt - theta_avg);
	}
	B = B*ngr_cnt/(nchain - 1.0);
	if(W > 0.0) {
	  R = ((1.0 - (1.0/ngr_cnt))*W + B/ngr_cnt)/W;
	  if(R > 0) {
	    R = sqrt(R);
	    if(R > 1.2) {test_pass = 1; break;}
	  }
	}
      }
    }
    if(test_pass) {
      /* This did not pass the Gelman & Rubin test, restart the chains for
         storage */
      nindx = 0;
      cur_chain_length = 0;
      NGR_pass = 0;
    } else {
      NGR_pass++;
    }
  }
  *Nlinks = cur_chain_length;

  if(NGR_pass < NGR_PASS_STOP_LIMIT) {
    fprintf(stderr,"#WARNING: The DEMC chain did not converge. Consider re-running with a higher maximum link threshhold (the current value is %d) using the --maxlinks option.\n", maxlinks);
  }
  
  for(i=0; i < nchain; i++) free(bvector[i]);
  free(bvector);
  free(chi2vector);
  free(dvector);
  free(gvector);
  free(mdata);
  free(errvector);
	
  return;

}

void PrintDEMCChain(FILE *outfile, _ModelData *m, double *l, int N) {
  int i;
  fprintf(outfile,"#logrho dlR dMbol log(M)  log(R)  log(Teff) Mbol BCV  BCJ  BCH  BCK  absVmag absJmag absHmag absKmag AV   distmod -2lnL\n");
  fprintf(outfile,"#[1]    [2] [3]   [4]     [5]     [6]       [7]  [8]  [9]  [10] [11] [12]    [13]    [14]    [15]    [16] [17]    [18]\n");
  for(i=0; i < N; i += 1) {
    fprintf(outfile,"%.17g", m[i].logrho);
    fprintf(outfile," %.17g", m[i].dlR);
    fprintf(outfile," %.17g", m[i].dMbol);
    fprintf(outfile," %.17g", m[i].lM);
    fprintf(outfile," %.17g", m[i].lR);
    fprintf(outfile," %.17g", m[i].lTeff);
    fprintf(outfile," %.17g", m[i].Mbol);
    fprintf(outfile," %.17g", m[i].BCV);
    fprintf(outfile," %.17g", m[i].BCJ);
    fprintf(outfile," %.17g", m[i].BCH);
    fprintf(outfile," %.17g", m[i].BCK);
    fprintf(outfile," %.17g", m[i].absVmag);
    fprintf(outfile," %.17g", m[i].absJmag);
    fprintf(outfile," %.17g", m[i].absHmag);
    fprintf(outfile," %.17g", m[i].absKmag);
    fprintf(outfile," %.17g", m[i].AV);
    fprintf(outfile," %.17g", m[i].distmod);
    fprintf(outfile," %.17g\n", l[i]);
  }
}

void PrintAutoFormat(char *s, double val, double err) {
  int edig;
  char F[256];
  char format[256];
  char v1s[256], v2s[256];
  int v1, v2, v1_, v2_;
  double v1d, v2d;
  if(err > 0) {
    edig = floor(log(err)/log(10.0)) - 1;
  } else {
    edig = 0;
  }
  if(edig < 0) {
    sprintf(F,".%df",-edig);
    sprintf(format,"%%s:\t %%%s +- %%%s\n", F, F);
    printf(format, s, val, err);
  }
  else {
    v1d = val*pow(10,-edig);
    v2d = err*pow(10,-edig);
    sprintf(v1s,"%.0f",v1d);
    sprintf(v2s,"%.0f",v2d);
    v1 = atoi(v1s);
    v2 = atoi(v2s);
    v1_ = v1*pow(10,edig);
    v2_ = v2*pow(10,edig);
    printf("%s:\t %d +- %d\n", s, v1_, v2_);
  }
}

void PrintOutputParameters(_ModelData *m,int N, double fracburnin, int isAV, int isdist)
{
  int i, istart;
  double *x, val1, val2;
  if((x = (double *) malloc(N * sizeof(double))) == NULL) {
    fprintf(stderr,"Memory Allocation Error\n");
    exit(1);
  }
  printf("#Parameters Determined From Empirical Model\n");
  printf("#Parameter:   MedianValue  +-   1_sigma_Error\n");

  istart = ceil(fracburnin*N);

  /* logrho */
  for(i=istart; i < N; i++){
    x[i-istart] = m[i].logrho;
  }
  val1 = ComputeMedian(x, N-istart);
  val2 = ComputeStddev(x, N-istart);
  PrintAutoFormat("log10(rho)", val1, val2);
  
  /* rho */
  for(i=istart; i < N; i++) {
    x[i-istart] = exp(M_LN10*m[i].logrho);
  }
  val1 = ComputeMedian(x, N-istart);
  val2 = ComputeStddev(x, N-istart);
  PrintAutoFormat("rho[cgs]", val1, val2);

  /* dlR */
  for(i=istart; i < N; i++) {
    x[i-istart] = m[i].dlR;
  }
  val1 = ComputeMedian(x, N-istart);
  val2 = ComputeStddev(x, N-istart);
  PrintAutoFormat("dlR", val1, val2);

  /* dMbol */
  for(i=istart; i < N; i++) {
    x[i-istart] = m[i].dMbol;
  }
  val1 = ComputeMedian(x, N-istart);
  val2 = ComputeStddev(x, N-istart);
  PrintAutoFormat("dMbol", val1, val2);

  /* lM */
  for(i=istart; i < N; i++) {
    x[i-istart] = m[i].lM;
  }
  val1 = ComputeMedian(x, N-istart);
  val2 = ComputeStddev(x, N-istart);
  PrintAutoFormat("log10(Mass)", val1, val2);

  /* Mass */
  for(i=istart; i < N; i++) {
    x[i-istart] = exp(M_LN10*m[i].lM);
  }
  val1 = ComputeMedian(x, N-istart);
  val2 = ComputeStddev(x, N-istart);
  PrintAutoFormat("Mass[Msun]", val1, val2);

  /* lR */
  for(i=istart; i < N; i++) {
    x[i-istart] = m[i].lR;
  }
  val1 = ComputeMedian(x, N-istart);
  val2 = ComputeStddev(x, N-istart);
  PrintAutoFormat("log10(Radius)", val1, val2);

  /* Radius */
  for(i=istart; i < N; i++) {
    x[i-istart] = exp(M_LN10*m[i].lR);
  }
  val1 = ComputeMedian(x, N-istart);
  val2 = ComputeStddev(x, N-istart);
  PrintAutoFormat("Radius[Rsun]", val1, val2);

  /* lTeff */
  for(i=istart; i < N; i++) {
    x[i-istart] = m[i].lTeff;
  }
  val1 = ComputeMedian(x, N-istart);
  val2 = ComputeStddev(x, N-istart);
  PrintAutoFormat("log10(Teff[K])", val1, val2);

  /* Teff */
  for(i=istart; i < N; i++) {
    x[i-istart] = exp(M_LN10*m[i].lTeff);
  }
  val1 = ComputeMedian(x, N-istart);
  val2 = ComputeStddev(x, N-istart);
  PrintAutoFormat("Teff[K]", val1, val2);

  /* Mbol */
  for(i=istart; i < N; i++) {
    x[i-istart] = m[i].Mbol;
  }
  val1 = ComputeMedian(x, N-istart);
  val2 = ComputeStddev(x, N-istart);
  PrintAutoFormat("Mbol", val1, val2);

  /* Abs MV */
  for(i=istart; i < N; i++) {
    x[i-istart] = m[i].absVmag;
  }
  val1 = ComputeMedian(x, N-istart);
  val2 = ComputeStddev(x, N-istart);
  PrintAutoFormat("absVmag", val1, val2);

  /* Abs MJ */
  for(i=istart; i < N; i++) {
    x[i-istart] = m[i].absJmag;
  }
  val1 = ComputeMedian(x, N-istart);
  val2 = ComputeStddev(x, N-istart);
  PrintAutoFormat("absJmag", val1, val2);

  /* Abs MH */
  for(i=istart; i < N; i++) {
    x[i-istart] = m[i].absHmag;
  }
  val1 = ComputeMedian(x, N-istart);
  val2 = ComputeStddev(x, N-istart);
  PrintAutoFormat("absHmag", val1, val2);

  /* Abs MK */
  for(i=istart; i < N; i++) {
    x[i-istart] = m[i].absKmag;
  }
  val1 = ComputeMedian(x, N-istart);
  val2 = ComputeStddev(x, N-istart);
  PrintAutoFormat("absKmag", val1, val2);

  if(isdist) {
    /* Apparent V */
    for(i=istart; i < N; i++) {
      x[i-istart] = m[i].absVmag + m[i].distmod + m[i].AV;
    }
    val1 = ComputeMedian(x, N-istart);
    val2 = ComputeStddev(x, N-istart);
    PrintAutoFormat("Vmag", val1, val2);
    
    /* Apparent J */
    for(i=istart; i < N; i++) {
      x[i-istart] = m[i].absJmag + m[i].distmod + m[i].AV*REDLAW_AJ;
    }
    val1 = ComputeMedian(x, N-istart);
    val2 = ComputeStddev(x, N-istart);
    PrintAutoFormat("Jmag", val1, val2);

    /* Apparent H */
    for(i=istart; i < N; i++) {
      x[i-istart] = m[i].absHmag + m[i].distmod + m[i].AV*REDLAW_AH;
    }
    val1 = ComputeMedian(x, N-istart);
    val2 = ComputeStddev(x, N-istart);
    PrintAutoFormat("Hmag", val1, val2);

    /* Apparent K */
    for(i=istart; i < N; i++) {
      x[i-istart] = m[i].absKmag + m[i].distmod + m[i].AV*REDLAW_AK;
    }
    val1 = ComputeMedian(x, N-istart);
    val2 = ComputeStddev(x, N-istart);
    PrintAutoFormat("Kmag", val1, val2);

    /* distance modulus */
    for(i=istart; i < N; i++) {
      x[i-istart] = m[i].distmod;
    }
    val1 = ComputeMedian(x, N-istart);
    val2 = ComputeStddev(x, N-istart);
    PrintAutoFormat("dist_modulus[mag]", val1, val2);

    /* distance */
    for(i=istart; i < N; i++) {
      x[i-istart] = exp(M_LN10*(m[i].distmod+5.)/5.);
    }
    val1 = ComputeMedian(x, N-istart);
    val2 = ComputeStddev(x, N-istart);
    PrintAutoFormat("distance[pc]", val1, val2);

    /* parallax */
    for(i=istart; i < N; i++) {
      x[i-istart] = 1000./exp(M_LN10*(m[i].distmod+5.)/5.);
    }
    val1 = ComputeMedian(x, N-istart);
    val2 = ComputeStddev(x, N-istart);
    PrintAutoFormat("parallax[mas]", val1, val2);

    /* Angular diameter */
    for(i=istart; i < N; i++) {
      x[i-istart] = 9.298438*exp(M_LN10*m[i].lR)/exp(M_LN10*((m[i].distmod + 5.0)/5.0));

    }
    val1 = ComputeMedian(x, N-istart);
    val2 = ComputeStddev(x, N-istart);
    PrintAutoFormat("Angular_Diameter[mas]", val1, val2);
  }

  /* Apparent V - J*/
  for(i=istart; i < N; i++) {
    x[i-istart] = m[i].absVmag - m[i].absJmag + m[i].AV*(1 - REDLAW_AJ);
  }
  val1 = ComputeMedian(x, N-istart);
  val2 = ComputeStddev(x, N-istart);
  PrintAutoFormat("V-J", val1, val2);

  /* Apparent V - H*/
  for(i=istart; i < N; i++) {
    x[i-istart] = m[i].absVmag - m[i].absHmag + m[i].AV*(1 - REDLAW_AH);
  }
  val1 = ComputeMedian(x, N-istart);
  val2 = ComputeStddev(x, N-istart);
  PrintAutoFormat("V-H", val1, val2);

  /* Apparent V - K*/
  for(i=istart; i < N; i++) {
    x[i-istart] = m[i].absVmag - m[i].absKmag + m[i].AV*(1 - REDLAW_AK);
  }
  val1 = ComputeMedian(x, N-istart);
  val2 = ComputeStddev(x, N-istart);
  PrintAutoFormat("V-K", val1, val2);

  /* Apparent J - H*/
  for(i=istart; i < N; i++) {
    x[i-istart] = m[i].absJmag - m[i].absHmag + m[i].AV*(REDLAW_AJ - REDLAW_AH);
  }
  val1 = ComputeMedian(x, N-istart);
  val2 = ComputeStddev(x, N-istart);
  PrintAutoFormat("J-H", val1, val2);

  /* Apparent J - K*/
  for(i=istart; i < N; i++) {
    x[i-istart] = m[i].absJmag - m[i].absKmag + m[i].AV*(REDLAW_AJ - REDLAW_AK);
  }
  val1 = ComputeMedian(x, N-istart);
  val2 = ComputeStddev(x, N-istart);
  PrintAutoFormat("J-K", val1, val2);

  /* Apparent H - K*/
  for(i=istart; i < N; i++) {
    x[i-istart] = m[i].absHmag - m[i].absKmag + m[i].AV*(REDLAW_AH - REDLAW_AK);
  }
  val1 = ComputeMedian(x, N-istart);
  val2 = ComputeStddev(x, N-istart);
  PrintAutoFormat("H-K", val1, val2);

  if(isAV) {
    /* AV */
    for(i=istart; i < N; i++) {
      x[i-istart] = m[i].AV;
    }
    val1 = ComputeMedian(x, N-istart);
    val2 = ComputeStddev(x, N-istart);
    PrintAutoFormat("AV[mag]", val1, val2);
  }
}

void Help(char *argv) {
  fprintf(stderr,
"\n"
"The code fits a set of observational constraints for a K and M dwarf star\n"
"using this empirical model. It runs a Differential Evolution Markov Chain\n"
"with the free parameters being log(rho), dlR, dMbol, and the parameters\n"
"describing the empirical model. Prior constraints based on the fit presented\n"
"in Hartman et al. 2014 are applied to these latter parameters. These\n"
"parameters are allowed to vary so that uncertainties in the empirical\n"
"model are propogated into the uncertainties on the resulting stellar\n"
"parameters."
"\n\n"
"One or more observational constraints are supplied on the command-line\n"
"using the syntax:\n\n"
"\t--parameter observed_value:1sigma_uncertainty[min_value:max_value]\n\n"
"Here observed_value and 1sigma_uncertainty are required by the program\n"
"[min_value:max_value] are optional additional constraints used to indicate\n"
"a minimum and/or maximum allowed value for the parameter. Use [:max_value]\n"
"or [min_value:] if there is only an upper or lower constraint.\n"
"\n"
"The follow observational constraints are allowed:\n\n"
"\t-m | --mass\t\t\t--- Mass of the star in solar units\n\n"
"\t-r | --radius\t\t\t--- Radius of the star in soloar units\n\n"
"\t--rho | --density\t\t--- Stellar density in grams per cm^3\n\n"
"\t--logrho | --lrho\t\t--- Log10(density) in grams per cm^3\n\n"
"\t--rhosolar | --densitysolar\t--- Stellar density in solar units\n\n"
"\t--logrhosolar | --lrhosolar\t--- Log10(density) in solar units\n\n"
"\t-logg | --logg\t\t\t--- log10(surface grav.) in cgs units\n\n"
"\t-loggsolar | --loggsolar\t--- log10(surface grav.) in solar units\n\n"
"\t-t | --teff\t\t\t--- Surface temperature in K\n\n"
"\t-V | --Vmag\t\t\t--- Apparent V magnitude of the star\n\n"
"\t-J | --Jmag\t\t\t--- Apparent J magnitude of the star\n\n"
"\t-H | --Hmag\t\t\t--- Apparent H magnitude of the star\n\n"
"\t-K | --Kmag\t\t\t--- Apparent K magnitude of the star\n\n"
"\t-MV | --absVmag\t\t\t--- Absolute V magnitude of the star\n\n"
"\t-MJ | --absJmag\t\t\t--- Absolute J magnitude of the star\n\n"
"\t-MH | --absHmag\t\t\t--- Absolute H magnitude of the star\n\n"
"\t-MK | --absKmag\t\t\t--- Absolute K magnitude of the star\n\n"
"\t-VJ | --V-J\t\t\t--- Apparent V-J color of the star\n\n"
"\t-VH | --V-H\t\t\t--- Apparent V-H color of the star\n\n"
"\t-VK | --V-K\t\t\t--- Apparent V-K color of the star\n\n"
"\t-JH | --J-H\t\t\t--- Apparent J-H color of the star\n\n"
"\t-JK | --J-K\t\t\t--- Apparent J-K color of the star\n\n"
"\t-HK | --H-K\t\t\t--- Apparent H-K color of the star\n\n"
"\t-dm | --distmod\t\t\t--- Distance modulus in magnitudes\n\n"
"\t-d | --distance\t\t\t--- Distance to the star in parsecs\n\n"
"\t--plx | --parallax\t\t--- Parallax to the star in mas\n\n"
"\t--angdiam\t\t\t--- Angular diameter of the star in mas\n\n"
"\t-AV | --AV\t\t\t--- V-band extinction. A Cardelli 1989,\n"
"\t\t\t\t\t    RV=3.1 extinction law is assumed.\n\n"
"\t-EBV | --EBV\t\t\t--- Excess B-V color of star. A\n"
"\t\t\t\t\t    Cardelli 1989, RV=3.1 extinction\n"
"\t\t\t\t\t    law is assumed.\n\n"
"Independent observations of the same parameter are allowed by repeating\n"
"the parameter on the command-line.\n"
"\n"
"The following other options may be given on the command-line:\n\n"
"\t--outchain file\t\t--- Output the DEMC chain to the specified\n"
"\t\t\t\t    file.\n\n"
"\t--ntest NlinksTest\t--- Calculate the Gelman & Rubin convergence\n"
"\t\t\t\t    criterian after NlinksTest links.\n"
"\t\t\t\t    The program will then continue\n"
"\t\t\t\t    running until the criterion is met.\n"
"\t\t\t\t    Default is 10000.\n\n"
"\t--maxlinks MaxNLinks\t--- Stop the DEMC after MaxNlinks regardless of\n"
"\t\t\t\t    whether or not the convergence criterion has been\n"
"\t\t\t\t    met. By default this is 1000000. If a negative\n"
"\t\t\t\t    value is given then there is no limit.\n\n"
"\t--maxmem MaxMem\t\t--- The maximum amount of memory to be allocated.\n"
"\t\t\t\t    If the chain exceeds this memory limit, new links\n"
"\t\t\t\t    will be written over the start of the chain.\n"
"\t\t\t\t    By default the limit is 2GB.\n\n"
"\t--burn-in fracburnin\t--- Skip the first fracburnin fraction\n"
"\t\t\t\t    of the links in calculating the median\n"
"\t\t\t\t    parameter values and uncertainties.\n"
"\t\t\t\t    Default is 0.1.\n\n"
"\t--useAV maxAV\t\t--- Allow for extinction in the model. A\n"
"\t\t\t\t    maximum optical extinction along the line\n"
"\t\t\t\t    of sight, in magnitudes, must be\n"
"\t\t\t\t    supplied. By default extinction\n"
"\t\t\t\t    is ignored in fitting the model.\n\n"
"\t--seed randseed\t\t--- Change the seed for the random number\n"
"\t\t\t\t    generator. randseed must be a\n"
"\t\t\t\t    (type long) integer.\n\n"
"\t--help\t\t\t--- Display this help.\n\n"
"Example:\n"
"--------\n\n"
"1. Estimate the parameters for GJ 1214\n"
"\t%s -V 14.64:0.04 -J 9.750:0.024 -H 9.094:0.024 \\\n"
"\t\t-K 8.782:0.020 -plx 68.71:0.6 --rho 27.7:8.9 \n\n"
"2. Estimate the parameters for a star with V-K = 3.990 +- 0.058\n"
"   rho = 3.88 +- 0.27 g/cm^3\n"
"\t%s --V-K 3.990:0.058 --rho 5.6:2.8\n\n"
	  , argv, argv);
  exit(1);
}

void Usage_noexit(char *argv) {
  fprintf(stderr,
	  "Usage: %s\n"
	  "\t<observed parameter list>\n"
	  "\t[--outchain file] [--ntest NlinksTest] [--maxlinks MaxNlinks]\n"
          "\t[--maxmem maxmem] [--burn-in fracburnin]\n"
	  "\t[--useAV maxAV] [--seed randseed] [--help]\n\n"
	  "This program estimates the physical properties of K and M dwarf stars\n"
	  "given an input set of observations, using the empirical stellar model\n"
	  "described in the paper Hartman et al. 2014, in preparation.\n"
, argv);
}

void Usage(char *argv) {
  Usage_noexit(argv);
  fprintf(stderr,
	  "Do \"%s --help\" for a list of observed parameters.\n\n", argv);
  exit(1);
}

int main(int argc, char **argv) {
  int i;
  _ModelData *OutChains;
  double *outloglike;
  int Nlinks;
  int useAV = 0;
  double maxAV = 0.;
  double fracburnin = 0.2;
  int paramcheck = 0;
  int nlinkstest = 10000;
  int maxlinks = 1000000;
  double maxmem = 2.0;
  int doprintchain = 0;
  int isdistmod, isAV;
  long seed = 0;

  FILE *outchainfile;

  _ObsData ObsData;
  ObsData.Nobs = 0;

  i=1;

  /* Parse the observed parameters from the command line */
  while(i < argc) {
    if(!strcmp(argv[i],"-m") || !strcmp(argv[i],"--mass") ||
       !strcmp(argv[i],"--m") || !strcmp(argv[i],"-mass")) {
      i++;
      if(i >= argc) Usage(argv[0]);
      parseobsdata(argv[0],argv[i-1],&ObsData,argv[i],OBS_TYPE_MASS);
      if(ObsData.obsval[ObsData.Nobs-1] > MAX_INPUT_MASS) {
	fprintf(stderr,"Error: this model may only be used for stars with masses less than 0.8 Msun.\n");
	exit(1);
      }
    }
    else if(!strcmp(argv[i],"-r") || !strcmp(argv[i],"--radius") ||
	    !strcmp(argv[i],"--r") || !strcmp(argv[i],"-radius")) {
      i++;
      if(i >= argc) Usage(argv[0]);
      parseobsdata(argv[0],argv[i-1],&ObsData,argv[i],OBS_TYPE_RADIUS);
      if(ObsData.obsval[ObsData.Nobs-1] > MAX_INPUT_RADIUS) {
	fprintf(stderr,"Error: this model may only be used for stars with radii less than 1.0 Rsun.\n");
	exit(1);
      }
    }
    else if(!strcmp(argv[i],"--rho") || !strcmp(argv[i],"--density") ||
	    !strcmp(argv[i],"-rho") || !strcmp(argv[i],"-density")) {
      i++;
      if(i >= argc) Usage(argv[0]);
      parseobsdata(argv[0],argv[i-1],&ObsData,argv[i],OBS_TYPE_DENSITY);
    }
    else if(!strcmp(argv[i],"--logrho") || !strcmp(argv[i],"--logdensity") ||
	    !strcmp(argv[i],"--lrho") || !strcmp(argv[i],"--ldensity")) {
      i++;
      if(i >= argc) Usage(argv[0]);
      parseobsdata(argv[0],argv[i-1],&ObsData,argv[i],OBS_TYPE_LOGDENSITY);
    }
    else if(!strcmp(argv[i],"-t") || !strcmp(argv[i],"--teff") ||
	    !strcmp(argv[i],"--t") || !strcmp(argv[i],"-teff")) {
      i++;
      if(i >= argc) Usage(argv[0]);
      parseobsdata(argv[0],argv[i-1],&ObsData,argv[i],OBS_TYPE_TEFF);
    }
    else if(!strcmp(argv[i],"-V") || !strcmp(argv[i],"-v") ||
	    !strcmp(argv[i],"--Vmag") || !strcmp(argv[i],"--vmag")) {
      i++;
      if(i >= argc) Usage(argv[0]);
      parseobsdata(argv[0],argv[i-1],&ObsData,argv[i],OBS_TYPE_VMAG);
    }
    else if(!strcmp(argv[i],"-J") || !strcmp(argv[i],"-j") ||
	    !strcmp(argv[i],"--Jmag") || !strcmp(argv[i],"--jmag")) {
      i++;
      if(i >= argc) Usage(argv[0]);
      parseobsdata(argv[0],argv[i-1],&ObsData,argv[i],OBS_TYPE_JMAG);
    }
    else if(!strcmp(argv[i],"-H") || !strcmp(argv[i],"-h") ||
	    !strcmp(argv[i],"--Hmag") || !strcmp(argv[i],"--hmag")) {
      i++;
      if(i >= argc) Usage(argv[0]);
      parseobsdata(argv[0],argv[i-1],&ObsData,argv[i],OBS_TYPE_HMAG);
    }
    else if(!strcmp(argv[i],"-K") || !strcmp(argv[i],"-k") ||
	    !strcmp(argv[i],"--Kmag") || !strcmp(argv[i],"--kmag")) {
      i++;
      if(i >= argc) Usage(argv[0]);
      parseobsdata(argv[0],argv[i-1],&ObsData,argv[i],OBS_TYPE_KMAG);
    }
    else if(!strcmp(argv[i],"-MV") || !strcmp(argv[i],"-Mv") ||
	    !strcmp(argv[i],"--MV") || !strcmp(argv[i],"--Mv") ||
	    !strcmp(argv[i],"--absVmag") || !strcmp(argv[i],"--absvmag")) {
      i++;
      if(i >= argc) Usage(argv[0]);
      parseobsdata(argv[0],argv[i-1],&ObsData,argv[i],OBS_TYPE_ABSVMAG);
    }
    else if(!strcmp(argv[i],"-MJ") || !strcmp(argv[i],"-Mj") ||
	    !strcmp(argv[i],"--MJ") || !strcmp(argv[i],"--Mj") ||
	    !strcmp(argv[i],"--absJmag") || !strcmp(argv[i],"--absjmag")) {
      i++;
      if(i >= argc) Usage(argv[0]);
      parseobsdata(argv[0],argv[i-1],&ObsData,argv[i],OBS_TYPE_ABSJMAG);
    }
    else if(!strcmp(argv[i],"-MH") || !strcmp(argv[i],"-Mh") ||
	    !strcmp(argv[i],"--MH") || !strcmp(argv[i],"--Mh") ||
	    !strcmp(argv[i],"--absHmag") || !strcmp(argv[i],"--abshmag")) {
      i++;
      if(i >= argc) Usage(argv[0]);
      parseobsdata(argv[0],argv[i-1],&ObsData,argv[i],OBS_TYPE_ABSHMAG);
    }
    else if(!strcmp(argv[i],"-MK") || !strcmp(argv[i],"-Mk") ||
	    !strcmp(argv[i],"--MK") || !strcmp(argv[i],"--Mk") ||
	    !strcmp(argv[i],"--absKmag") || !strcmp(argv[i],"--abskmag")) {
      i++;
      if(i >= argc) Usage(argv[0]);
      parseobsdata(argv[0],argv[i-1],&ObsData,argv[i],OBS_TYPE_ABSKMAG);
    }
    else if(!strcmp(argv[i],"-VJ") || !strcmp(argv[i],"-vj") ||
	    !strcmp(argv[i],"--VJ") || !strcmp(argv[i],"--vj") ||
	    !strcmp(argv[i],"--V-J") || !strcmp(argv[i],"--v-j") ||
	    !strcmp(argv[i],"-V-J") || !strcmp(argv[i],"-v-j")) {
      i++;
      if(i >= argc) Usage(argv[0]);
      parseobsdata(argv[0],argv[i-1],&ObsData,argv[i],OBS_TYPE_VJCOLOR);
    }
    else if(!strcmp(argv[i],"-VH") || !strcmp(argv[i],"-vh") ||
	    !strcmp(argv[i],"--VH") || !strcmp(argv[i],"--vh") ||
	    !strcmp(argv[i],"--V-H") || !strcmp(argv[i],"--v-h") ||
	    !strcmp(argv[i],"-V-H") || !strcmp(argv[i],"-v-h")) {
      i++;
      if(i >= argc) Usage(argv[0]);
      parseobsdata(argv[0],argv[i-1],&ObsData,argv[i],OBS_TYPE_VHCOLOR);
    }
    else if(!strcmp(argv[i],"-VK") || !strcmp(argv[i],"-vk") ||
	    !strcmp(argv[i],"--VK") || !strcmp(argv[i],"--vk") ||
	    !strcmp(argv[i],"--V-K") || !strcmp(argv[i],"--v-k") ||
	    !strcmp(argv[i],"-V-K") || !strcmp(argv[i],"-v-k")) {
      i++;
      if(i >= argc) Usage(argv[0]);
      parseobsdata(argv[0],argv[i-1],&ObsData,argv[i],OBS_TYPE_VKCOLOR);
    }
    else if(!strcmp(argv[i],"-JH") || !strcmp(argv[i],"-jh") ||
	    !strcmp(argv[i],"--JH") || !strcmp(argv[i],"--jh") ||
	    !strcmp(argv[i],"--J-H") || !strcmp(argv[i],"--j-h") ||
	    !strcmp(argv[i],"-J-H") || !strcmp(argv[i],"-j-h")) {
      i++;
      if(i >= argc) Usage(argv[0]);
      parseobsdata(argv[0],argv[i-1],&ObsData,argv[i],OBS_TYPE_JHCOLOR);
    }
    else if(!strcmp(argv[i],"-JK") || !strcmp(argv[i],"-jk") ||
	    !strcmp(argv[i],"--JK") || !strcmp(argv[i],"--jk") ||
	    !strcmp(argv[i],"--J-K") || !strcmp(argv[i],"--j-k") ||
	    !strcmp(argv[i],"-J-K") || !strcmp(argv[i],"-j-k")) {
      i++;
      if(i >= argc) Usage(argv[0]);
      parseobsdata(argv[0],argv[i-1],&ObsData,argv[i],OBS_TYPE_JKCOLOR);
    }
    else if(!strcmp(argv[i],"-HK") || !strcmp(argv[i],"-hk") ||
	    !strcmp(argv[i],"--HK") || !strcmp(argv[i],"--hk") ||
	    !strcmp(argv[i],"--H-K") || !strcmp(argv[i],"--h-k") ||
	    !strcmp(argv[i],"-H-K") || !strcmp(argv[i],"-h-k")) {
      i++;
      if(i >= argc) Usage(argv[0]);
      parseobsdata(argv[0],argv[i-1],&ObsData,argv[i],OBS_TYPE_HKCOLOR);
    }
    else if(!strcmp(argv[i],"-DM") || !strcmp(argv[i],"-dm") ||
	    !strcmp(argv[i],"--DM") || !strcmp(argv[i],"--dm") ||
	    !strcmp(argv[i],"--distmod") || !strcmp(argv[i],"--Distmod")) {
      i++;
      if(i >= argc) Usage(argv[0]);
      parseobsdata(argv[0],argv[i-1],&ObsData,argv[i],OBS_TYPE_DISTMOD);
    }
    else if(!strcmp(argv[i],"-D") || !strcmp(argv[i],"-d") ||
	    !strcmp(argv[i],"--Dist") || !strcmp(argv[i],"--dist") ||
	    !strcmp(argv[i],"--distance")) {
      i++;
      if(i >= argc) Usage(argv[0]);
      parseobsdata(argv[0],argv[i-1],&ObsData,argv[i],OBS_TYPE_DIST);
    }
    else if(!strcmp(argv[i],"-plx") || !strcmp(argv[i],"--plx") ||
	    !strcmp(argv[i],"--parallax") || !strcmp(argv[i],"-parallax")){
      i++;
      if(i >= argc) Usage(argv[0]);
      parseobsdata(argv[0],argv[i-1],&ObsData,argv[i],OBS_TYPE_PARALLAX);
    }
    else if(!strcmp(argv[i],"--angdiam")){
      i++;
      if(i >= argc) Usage(argv[0]);
      parseobsdata(argv[0],argv[i-1],&ObsData,argv[i],OBS_TYPE_ANGDIAM);
    }
    else if(!strcmp(argv[i],"-AV") || !strcmp(argv[i],"-av") ||
	    !strcmp(argv[i],"--AV") || !strcmp(argv[i],"--av")) {
      i++;
      if(i >= argc) Usage(argv[0]);
      parseobsdata(argv[0],argv[i-1],&ObsData,argv[i],OBS_TYPE_AV);
    }
    else if(!strcmp(argv[i],"-EBV") || !strcmp(argv[i],"-ebv") ||
	    !strcmp(argv[i],"--EBV") || !strcmp(argv[i],"--ebv") ||
	    !strcmp(argv[i],"-E(B-V)") || !strcmp(argv[i],"--E(B-V)")) {
      i++;
      if(i >= argc) Usage(argv[0]);
      parseobsdata(argv[0],argv[i-1],&ObsData,argv[i],OBS_TYPE_EBV);
    }
    else if(!strcmp(argv[i],"-logg") || !strcmp(argv[i],"--logg")) {
      i++;
      if(i >= argc) Usage(argv[0]);
      parseobsdata(argv[0],argv[i-1],&ObsData,argv[i],OBS_TYPE_LOGG);
    }
    else if(!strcmp(argv[i],"-loggsolar") || !strcmp(argv[i],"--loggsolar")) {
      i++;
      if(i >= argc) Usage(argv[0]);
      parseobsdata(argv[0],argv[i-1],&ObsData,argv[i],OBS_TYPE_LOGGSOLAR);
    }
    else if(!strcmp(argv[i],"--rhosolar") || !strcmp(argv[i],"--densitysolar") ||
	    !strcmp(argv[i],"-rhosolar") || !strcmp(argv[i],"-densitysolar")) {
      i++;
      if(i >= argc) Usage(argv[0]);
      parseobsdata(argv[0],argv[i-1],&ObsData,argv[i],OBS_TYPE_DENSITYSOLAR);
    }
    else if(!strcmp(argv[i],"--logrhosolar") || !strcmp(argv[i],"--logdensitysolar") ||
	    !strcmp(argv[i],"--lrhosolar") || !strcmp(argv[i],"--ldensitysolar")) {
      i++;
      if(i >= argc) Usage(argv[0]);
      parseobsdata(argv[0],argv[i-1],&ObsData,argv[i],OBS_TYPE_LOGDENSITYSOLAR);
    }
    else if(!strcmp(argv[i],"--outchain")) {
      i++;
      doprintchain = 1;
      if(i >= argc) Usage(argv[0]);
      if((outchainfile = fopen(argv[i],"w")) == NULL) {
	fprintf(stderr,"Cannot Write To %s\n", argv[i]);
	exit(2);
      }
    }
    else if(!strcmp(argv[i],"--useAV")) {
      i++;
      useAV = 1;
      if(i >= argc) Usage(argv[0]);
      maxAV = atof(argv[i]);
    }
    else if(!strcmp(argv[i],"--ntest")) {
      i++;
      if(i >= argc) Usage(argv[0]);
      nlinkstest = atoi(argv[i]);
    }
    else if(!strcmp(argv[i],"--maxlinks")) {
      i++;
      if(i >= argc) Usage(argv[0]);
      maxlinks = atoi(argv[i]);
    }
    else if(!strcmp(argv[i],"--maxmem")) {
      i++;
      if(i >= argc) Usage(argv[0]);
      maxmem = atof(argv[i]);
    }
    else if(!strcmp(argv[i],"--burn-in")) {
      i++;
      if(i >= argc) Usage(argv[0]);
      fracburnin = atof(argv[i]);
    }
    else if(!strcmp(argv[i],"--seed")) {
      i++;
      if(i >= argc) Usage(argv[0]);
      seed = atol(argv[i]);
    }
    else if(!strcmp(argv[i],"--help")) {
      Usage_noexit(argv[0]);
      Help(argv[0]);
    }
    else {
      Usage(argv[0]);
    }
    i++;
  }

  if(ObsData.Nobs == 0) {
    Usage_noexit(argv[0]);
    fprintf(stderr,"\nNo input parameters given. Do \"%s --help\" for help.\n", argv[0]);
    exit(1);
  }

  /* Give a warning if the only parameters are NIR photometric values */
  paramcheck = 0;
  for(i=0; i < ObsData.Nobs; i++) {
    if(ObsData.obstype[i] != OBS_TYPE_JMAG &&
       ObsData.obstype[i] != OBS_TYPE_HMAG &&
       ObsData.obstype[i] != OBS_TYPE_KMAG &&
       ObsData.obstype[i] != OBS_TYPE_JHCOLOR &&
       ObsData.obstype[i] != OBS_TYPE_JKCOLOR &&
       ObsData.obstype[i] != OBS_TYPE_HKCOLOR &&
       ObsData.obstype[i] != OBS_TYPE_AV &&
       ObsData.obstype[i] != OBS_TYPE_EBV) {
      paramcheck = 1;
    }
  }

  if(!paramcheck) {
    fprintf(stderr,"#WARNING: the stellar parameters are unreliable when apparent NIR photometric measurements are the only observables used.\n");
  }

  srand48(seed);

  RunDEMCMC(&ObsData, &OutChains, &outloglike, &Nlinks, useAV, maxAV, nlinkstest,maxlinks,maxmem,&isdistmod,&isAV);
  if(doprintchain){
    PrintDEMCChain(outchainfile,OutChains,outloglike,Nlinks);
    fclose(outchainfile);
  }
  PrintOutputParameters(OutChains,Nlinks,fracburnin, isAV, isdistmod);
  exit(0);
}