#include "i3_model.h"
#include "i3_options.h"
#include "i3_model_tools.h"
#include "string.h"
#include <stdlib.h>
#include <stdio.h>


void i3_model_setup_offsets(i3_model * model){
	model->nbytes = 0;
	model->byte_offsets = malloc(sizeof(int)*model->nparam);
	if (!strcmp(model->name,"sersics")){
		i3_sersics_parameter_set P;
		model->nbytes = sizeof(P);
		model->byte_offsets[0] = ((off_type)&(P.x0)) - (off_type)&P;
		model->byte_offsets[1] = ((off_type)&(P.y0)) - (off_type)&P;
		model->byte_offsets[2] = ((off_type)&(P.e1)) - (off_type)&P;
		model->byte_offsets[3] = ((off_type)&(P.e2)) - (off_type)&P;
		model->byte_offsets[4] = ((off_type)&(P.radius)) - (off_type)&P;
		model->byte_offsets[5] = ((off_type)&(P.radius_ratio)) - (off_type)&P;
		model->byte_offsets[6] = ((off_type)&(P.bulge_A)) - (off_type)&P;
		model->byte_offsets[7] = ((off_type)&(P.disc_A)) - (off_type)&P;
		model->byte_offsets[8] = ((off_type)&(P.bulge_index)) - (off_type)&P;
		model->byte_offsets[9] = ((off_type)&(P.disc_index)) - (off_type)&P;
		model->byte_offsets[10] = ((off_type)&(P.delta_e_bulge)) - (off_type)&P;
		model->byte_offsets[11] = ((off_type)&(P.delta_theta_bulge)) - (off_type)&P;

}
	else if (!strcmp(model->name,"logsersics")){
		i3_logsersics_parameter_set P;
		model->nbytes = sizeof(P);
		model->byte_offsets[0] = ((off_type)&(P.x0)) - (off_type)&P;
		model->byte_offsets[1] = ((off_type)&(P.y0)) - (off_type)&P;
		model->byte_offsets[2] = ((off_type)&(P.e1)) - (off_type)&P;
		model->byte_offsets[3] = ((off_type)&(P.e2)) - (off_type)&P;
		model->byte_offsets[4] = ((off_type)&(P.radius)) - (off_type)&P;
		model->byte_offsets[5] = ((off_type)&(P.radius_ratio)) - (off_type)&P;
		model->byte_offsets[6] = ((off_type)&(P.bulge_A)) - (off_type)&P;
		model->byte_offsets[7] = ((off_type)&(P.disc_A)) - (off_type)&P;
		model->byte_offsets[8] = ((off_type)&(P.bulge_index)) - (off_type)&P;
		model->byte_offsets[9] = ((off_type)&(P.disc_index)) - (off_type)&P;
		model->byte_offsets[10] = ((off_type)&(P.delta_e_bulge)) - (off_type)&P;
		model->byte_offsets[11] = ((off_type)&(P.delta_theta_bulge)) - (off_type)&P;

}
	else if (!strcmp(model->name,"quadratic_test")){
		i3_quadratic_test_parameter_set P;
		model->nbytes = sizeof(P);
		model->byte_offsets[0] = ((off_type)&(P.x1)) - (off_type)&P;
		model->byte_offsets[1] = ((off_type)&(P.x2)) - (off_type)&P;
		model->byte_offsets[2] = ((off_type)&(P.x3)) - (off_type)&P;
		model->byte_offsets[3] = ((off_type)&(P.x4)) - (off_type)&P;
		model->byte_offsets[4] = ((off_type)&(P.x5)) - (off_type)&P;

}
	else if (!strcmp(model->name,"multisersics")){
		i3_multisersics_parameter_set P;
		model->nbytes = sizeof(P);
		model->byte_offsets[0] = ((off_type)&(P.ra_as)) - (off_type)&P;
		model->byte_offsets[1] = ((off_type)&(P.dec_as)) - (off_type)&P;
		model->byte_offsets[2] = ((off_type)&(P.e1)) - (off_type)&P;
		model->byte_offsets[3] = ((off_type)&(P.e2)) - (off_type)&P;
		model->byte_offsets[4] = ((off_type)&(P.radius)) - (off_type)&P;
		model->byte_offsets[5] = ((off_type)&(P.radius_ratio)) - (off_type)&P;
		model->byte_offsets[6] = ((off_type)&(P.bulge_A)) - (off_type)&P;
		model->byte_offsets[7] = ((off_type)&(P.disc_A)) - (off_type)&P;
		model->byte_offsets[8] = ((off_type)&(P.bulge_index)) - (off_type)&P;
		model->byte_offsets[9] = ((off_type)&(P.disc_index)) - (off_type)&P;

}
	else if (!strcmp(model->name,"multixray")){
		i3_multixray_parameter_set P;
		model->nbytes = sizeof(P);
		model->byte_offsets[0] = ((off_type)&(P.alpha_arcsec)) - (off_type)&P;
		model->byte_offsets[1] = ((off_type)&(P.delta_arcsec)) - (off_type)&P;
		model->byte_offsets[2] = ((off_type)&(P.e1)) - (off_type)&P;
		model->byte_offsets[3] = ((off_type)&(P.e2)) - (off_type)&P;
		model->byte_offsets[4] = ((off_type)&(P.rc)) - (off_type)&P;
		model->byte_offsets[5] = ((off_type)&(P.A)) - (off_type)&P;
		model->byte_offsets[6] = ((off_type)&(P.beta)) - (off_type)&P;
		model->byte_offsets[7] = ((off_type)&(P.background)) - (off_type)&P;

}
}

void i3_model_destroy(i3_model * model){
	free(model->param_fixed);
	free(model->byte_offsets);
	free(model->min);
	free(model->max);
	free(model->param_type);
	free(model);
}

i3_model * i3_model_create(char * name, i3_options * options){
	i3_model * model = malloc(sizeof(i3_model));
	strcpy(model->name,name);
	
	if (!strcmp(name,"sersics")) {
		model->nparam=i3_sersics_nparam;
		model->likelihood = (likelihood_function) &i3_sersics_likelihood;
		model->prior = (prior_function) NULL;
		model->propose = (proposal_function) &i3_covariance_matrix_proposal;
		model->posterior_derivative = (derivative_function) NULL;
		model->start = (start_position_function) &i3_sersics_start;
		model->map_physical = (map_physical_function) &i3_sersics_beermat_mapping;
		/* Set up array describing whether parameters are fixed. */
		model->param_fixed = malloc(sizeof(bool)*model->nparam);
		model->param_fixed[0] = options->sersics_x0_fixed;  /* The x-coordinate of the center of the de Vaucouleurs bulge profile*/
		model->param_fixed[1] = options->sersics_y0_fixed;  /* The y-coordinate of the center of the de Vaucouleurs bulge profile*/
		model->param_fixed[2] = options->sersics_e1_fixed;  /* The ellipticity along x axs*/
		model->param_fixed[3] = options->sersics_e2_fixed;  /* The ellipticity along x=y axs*/
		model->param_fixed[4] = options->sersics_radius_fixed;  /* The bulge scale radius (exp scale radius is set from this too)*/
		model->param_fixed[5] = options->sersics_radius_ratio_fixed;  /* The bulge-disc ratio*/
		model->param_fixed[6] = options->sersics_bulge_A_fixed;  /* amplitude of the bulge*/
		model->param_fixed[7] = options->sersics_disc_A_fixed;  /* amplitude of the disc*/
		model->param_fixed[8] = options->sersics_bulge_index_fixed;  /* Sersic index of the bulge component*/
		model->param_fixed[9] = options->sersics_disc_index_fixed;  /* Sersic index of the disc component*/
		model->param_fixed[10] = options->sersics_delta_e_bulge_fixed;  /* The additional e in the bulge component compared to the disc*/
		model->param_fixed[11] = options->sersics_delta_theta_bulge_fixed;  /* The rotation of the the bulge compared to the disc*/
		/*  Set up array of parameter types */
		model->param_type = malloc(sizeof(i3_parameter_type)*model->nparam);
		model->param_type[0] = i3_parameter_type_flt;
		model->param_type[1] = i3_parameter_type_flt;
		model->param_type[2] = i3_parameter_type_flt;
		model->param_type[3] = i3_parameter_type_flt;
		model->param_type[4] = i3_parameter_type_flt;
		model->param_type[5] = i3_parameter_type_flt;
		model->param_type[6] = i3_parameter_type_flt;
		model->param_type[7] = i3_parameter_type_flt;
		model->param_type[8] = i3_parameter_type_flt;
		model->param_type[9] = i3_parameter_type_flt;
		model->param_type[10] = i3_parameter_type_flt;
		model->param_type[11] = i3_parameter_type_flt;
		
		
	}
	
	else 
	if (!strcmp(name,"logsersics")) {
		model->nparam=i3_logsersics_nparam;
		model->likelihood = (likelihood_function) &i3_logsersics_likelihood;
		model->prior = (prior_function) NULL;
		model->propose = (proposal_function) &i3_covariance_matrix_proposal;
		model->posterior_derivative = (derivative_function) NULL;
		model->start = (start_position_function) &i3_logsersics_start;
		model->map_physical = (map_physical_function) &i3_logsersics_beermat_mapping;
		/* Set up array describing whether parameters are fixed. */
		model->param_fixed = malloc(sizeof(bool)*model->nparam);
		model->param_fixed[0] = options->logsersics_x0_fixed;  /* The x-coordinate of the center of the de Vaucouleurs bulge profile*/
		model->param_fixed[1] = options->logsersics_y0_fixed;  /* The y-coordinate of the center of the de Vaucouleurs bulge profile*/
		model->param_fixed[2] = options->logsersics_e1_fixed;  /* The ellipticity along x axs*/
		model->param_fixed[3] = options->logsersics_e2_fixed;  /* The ellipticity along x=y axs*/
		model->param_fixed[4] = options->logsersics_radius_fixed;  /* The bulge scale radius (exp scale radius is set from this too)*/
		model->param_fixed[5] = options->logsersics_radius_ratio_fixed;  /* The bulge-disc ratio*/
		model->param_fixed[6] = options->logsersics_bulge_A_fixed;  /* Log of the flux in the bulge*/
		model->param_fixed[7] = options->logsersics_disc_A_fixed;  /* Log of the flux in the disc*/
		model->param_fixed[8] = options->logsersics_bulge_index_fixed;  /* Sersic index of the bulge component*/
		model->param_fixed[9] = options->logsersics_disc_index_fixed;  /* Sersic index of the disc component*/
		model->param_fixed[10] = options->logsersics_delta_e_bulge_fixed;  /* The additional e in the bulge component compared to the disc*/
		model->param_fixed[11] = options->logsersics_delta_theta_bulge_fixed;  /* The rotation of the the bulge compared to the disc*/
		/*  Set up array of parameter types */
		model->param_type = malloc(sizeof(i3_parameter_type)*model->nparam);
		model->param_type[0] = i3_parameter_type_flt;
		model->param_type[1] = i3_parameter_type_flt;
		model->param_type[2] = i3_parameter_type_flt;
		model->param_type[3] = i3_parameter_type_flt;
		model->param_type[4] = i3_parameter_type_flt;
		model->param_type[5] = i3_parameter_type_flt;
		model->param_type[6] = i3_parameter_type_flt;
		model->param_type[7] = i3_parameter_type_flt;
		model->param_type[8] = i3_parameter_type_flt;
		model->param_type[9] = i3_parameter_type_flt;
		model->param_type[10] = i3_parameter_type_flt;
		model->param_type[11] = i3_parameter_type_flt;
		
		
	}
	
	else 
	if (!strcmp(name,"quadratic_test")) {
		model->nparam=i3_quadratic_test_nparam;
		model->likelihood = (likelihood_function) &i3_quadratic_test_likelihood;
		model->prior = (prior_function) NULL;
		model->propose = (proposal_function) &i3_covariance_matrix_proposal;
		model->posterior_derivative = (derivative_function) NULL;
		model->start = (start_position_function) NULL;
		model->map_physical = (map_physical_function) NULL;
		/* Set up array describing whether parameters are fixed. */
		model->param_fixed = malloc(sizeof(bool)*model->nparam);
		model->param_fixed[0] = options->quadratic_test_x1_fixed;  /* Center of gaussian in dimension 1*/
		model->param_fixed[1] = options->quadratic_test_x2_fixed;  /* Center of gaussian  in dimension 2*/
		model->param_fixed[2] = options->quadratic_test_x3_fixed;  /* Center of gaussian  in dimension 3*/
		model->param_fixed[3] = options->quadratic_test_x4_fixed;  /* Center of gaussian  in dimension 4*/
		model->param_fixed[4] = options->quadratic_test_x5_fixed;  /* Unused test float parameter*/
		/*  Set up array of parameter types */
		model->param_type = malloc(sizeof(i3_parameter_type)*model->nparam);
		model->param_type[0] = i3_parameter_type_flt;
		model->param_type[1] = i3_parameter_type_flt;
		model->param_type[2] = i3_parameter_type_flt;
		model->param_type[3] = i3_parameter_type_flt;
		model->param_type[4] = i3_parameter_type_array;
		
		
	}
	
	else 
	if (!strcmp(name,"multisersics")) {
		model->nparam=i3_multisersics_nparam;
		model->likelihood = (likelihood_function) &i3_multisersics_likelihood;
		model->prior = (prior_function) NULL;
		model->propose = (proposal_function) &i3_covariance_matrix_proposal;
		model->posterior_derivative = (derivative_function) NULL;
		model->start = (start_position_function) &i3_multisersics_start;
		model->map_physical = (map_physical_function) &i3_multisersics_beermat_mapping;
		/* Set up array describing whether parameters are fixed. */
		model->param_fixed = malloc(sizeof(bool)*model->nparam);
		model->param_fixed[0] = options->multisersics_ra_as_fixed;  /* The x-coordinate of the center of the de Vaucouleurs bulge profile*/
		model->param_fixed[1] = options->multisersics_dec_as_fixed;  /* The y-coordinate of the center of the de Vaucouleurs bulge profile*/
		model->param_fixed[2] = options->multisersics_e1_fixed;  /* The ellipticity along x axs*/
		model->param_fixed[3] = options->multisersics_e2_fixed;  /* The ellipticity along x=y axs*/
		model->param_fixed[4] = options->multisersics_radius_fixed;  /* The bulge scale radius (exp scale radius is set from this too)*/
		model->param_fixed[5] = options->multisersics_radius_ratio_fixed;  /* The bulge-disc ratio*/
		model->param_fixed[6] = options->multisersics_bulge_A_fixed;  /* amplitude of the bulge*/
		model->param_fixed[7] = options->multisersics_disc_A_fixed;  /* amplitude of the disc*/
		model->param_fixed[8] = options->multisersics_bulge_index_fixed;  /* Sersic index of the bulge component*/
		model->param_fixed[9] = options->multisersics_disc_index_fixed;  /* Sersic index of the disc component*/
		/*  Set up array of parameter types */
		model->param_type = malloc(sizeof(i3_parameter_type)*model->nparam);
		model->param_type[0] = i3_parameter_type_flt;
		model->param_type[1] = i3_parameter_type_flt;
		model->param_type[2] = i3_parameter_type_flt;
		model->param_type[3] = i3_parameter_type_flt;
		model->param_type[4] = i3_parameter_type_flt;
		model->param_type[5] = i3_parameter_type_flt;
		model->param_type[6] = i3_parameter_type_flt;
		model->param_type[7] = i3_parameter_type_flt;
		model->param_type[8] = i3_parameter_type_flt;
		model->param_type[9] = i3_parameter_type_flt;
		
		
	}
	
	else 
	if (!strcmp(name,"multixray")) {
		model->nparam=i3_multixray_nparam;
		model->likelihood = (likelihood_function) &i3_multixray_likelihood;
		model->prior = (prior_function) NULL;
		model->propose = (proposal_function) &i3_covariance_matrix_proposal;
		model->posterior_derivative = (derivative_function) NULL;
		model->start = (start_position_function) NULL;
		model->map_physical = (map_physical_function) &i3_multixray_map;
		/* Set up array describing whether parameters are fixed. */
		model->param_fixed = malloc(sizeof(bool)*model->nparam);
		model->param_fixed[0] = options->multixray_alpha_arcsec_fixed;  /* Tangent plane RA offset of centroid from baseline*/
		model->param_fixed[1] = options->multixray_delta_arcsec_fixed;  /* Tangent plane DEC offset of centroid from baseline*/
		model->param_fixed[2] = options->multixray_e1_fixed;  /* The ellipticity along x axs*/
		model->param_fixed[3] = options->multixray_e2_fixed;  /* The ellipticity along x=y axs*/
		model->param_fixed[4] = options->multixray_rc_fixed;  /* The outer radius r_c*/
		model->param_fixed[5] = options->multixray_A_fixed;  /* Amplitude of the profile*/
		model->param_fixed[6] = options->multixray_beta_fixed;  
		model->param_fixed[7] = options->multixray_background_fixed;  /* A background noise level to add to the model*/
		/*  Set up array of parameter types */
		model->param_type = malloc(sizeof(i3_parameter_type)*model->nparam);
		model->param_type[0] = i3_parameter_type_flt;
		model->param_type[1] = i3_parameter_type_flt;
		model->param_type[2] = i3_parameter_type_flt;
		model->param_type[3] = i3_parameter_type_flt;
		model->param_type[4] = i3_parameter_type_flt;
		model->param_type[5] = i3_parameter_type_flt;
		model->param_type[6] = i3_parameter_type_flt;
		model->param_type[7] = i3_parameter_type_array;
		
		
	}
	
	else{
		free(model);
		return NULL;
	}
	i3_model_setup_offsets(model);
	model->min = i3_model_option_min(model,options);
	model->max = i3_model_option_max(model,options);
	return model;	
}

i3_parameter_set * i3_model_option_widths(i3_model * model, i3_options * options){
	i3_parameter_set * p = (i3_parameter_set*)malloc(model->nbytes);
	if (!strcmp(model->name,"sersics")){
		i3_sersics_parameter_set * P = (i3_sersics_parameter_set*) p;
		P->x0=options->sersics_x0_width;
		P->y0=options->sersics_y0_width;
		P->e1=options->sersics_e1_width;
		P->e2=options->sersics_e2_width;
		P->radius=options->sersics_radius_width;
		P->radius_ratio=options->sersics_radius_ratio_width;
		P->bulge_A=options->sersics_bulge_A_width;
		P->disc_A=options->sersics_disc_A_width;
		P->bulge_index=options->sersics_bulge_index_width;
		P->disc_index=options->sersics_disc_index_width;
		P->delta_e_bulge=options->sersics_delta_e_bulge_width;
		P->delta_theta_bulge=options->sersics_delta_theta_bulge_width;
}
	else if (!strcmp(model->name,"logsersics")){
		i3_logsersics_parameter_set * P = (i3_logsersics_parameter_set*) p;
		P->x0=options->logsersics_x0_width;
		P->y0=options->logsersics_y0_width;
		P->e1=options->logsersics_e1_width;
		P->e2=options->logsersics_e2_width;
		P->radius=options->logsersics_radius_width;
		P->radius_ratio=options->logsersics_radius_ratio_width;
		P->bulge_A=options->logsersics_bulge_A_width;
		P->disc_A=options->logsersics_disc_A_width;
		P->bulge_index=options->logsersics_bulge_index_width;
		P->disc_index=options->logsersics_disc_index_width;
		P->delta_e_bulge=options->logsersics_delta_e_bulge_width;
		P->delta_theta_bulge=options->logsersics_delta_theta_bulge_width;
}
	else if (!strcmp(model->name,"quadratic_test")){
		i3_quadratic_test_parameter_set * P = (i3_quadratic_test_parameter_set*) p;
		P->x1=options->quadratic_test_x1_width;
		P->x2=options->quadratic_test_x2_width;
		P->x3=options->quadratic_test_x3_width;
		P->x4=options->quadratic_test_x4_width;
		for (int ip=0; ip<3; ip++) P->x5[ip]=options->quadratic_test_x5_width[ip];
}
	else if (!strcmp(model->name,"multisersics")){
		i3_multisersics_parameter_set * P = (i3_multisersics_parameter_set*) p;
		P->ra_as=options->multisersics_ra_as_width;
		P->dec_as=options->multisersics_dec_as_width;
		P->e1=options->multisersics_e1_width;
		P->e2=options->multisersics_e2_width;
		P->radius=options->multisersics_radius_width;
		P->radius_ratio=options->multisersics_radius_ratio_width;
		P->bulge_A=options->multisersics_bulge_A_width;
		P->disc_A=options->multisersics_disc_A_width;
		P->bulge_index=options->multisersics_bulge_index_width;
		P->disc_index=options->multisersics_disc_index_width;
}
	else if (!strcmp(model->name,"multixray")){
		i3_multixray_parameter_set * P = (i3_multixray_parameter_set*) p;
		P->alpha_arcsec=options->multixray_alpha_arcsec_width;
		P->delta_arcsec=options->multixray_delta_arcsec_width;
		P->e1=options->multixray_e1_width;
		P->e2=options->multixray_e2_width;
		P->rc=options->multixray_rc_width;
		P->A=options->multixray_A_width;
		P->beta=options->multixray_beta_width;
		for (int ip=0; ip<30; ip++) P->background[ip]=options->multixray_background_width[ip];
}
	return p;
}

void i3_model_extract_ellipticity(i3_model * model, i3_parameter_set * p, i3_flt * e1, i3_flt * e2){
		if (!strcmp(model->name,"sersics")){
		i3_sersics_parameter_set * P = (i3_sersics_parameter_set*) p;
		 * e1 = P->e1;
		 * e2 = P->e1;
	}
	if (!strcmp(model->name,"logsersics")){
		i3_logsersics_parameter_set * P = (i3_logsersics_parameter_set*) p;
		 * e1 = P->e1;
		 * e2 = P->e1;
	}
	if (!strcmp(model->name,"quadratic_test")){
		I3_FATAL("Tried to extract e1 and e2 from parameter set for model 'quadratic_test', which has no e1 or e2 component",1);
	}
	if (!strcmp(model->name,"multisersics")){
		i3_multisersics_parameter_set * P = (i3_multisersics_parameter_set*) p;
		 * e1 = P->e1;
		 * e2 = P->e1;
	}
	if (!strcmp(model->name,"multixray")){
		i3_multixray_parameter_set * P = (i3_multixray_parameter_set*) p;
		 * e1 = P->e1;
		 * e2 = P->e1;
	}
}

bool i3_model_any_nan(i3_model * model, i3_parameter_set * p)
{
		if (!strcmp(model->name,"sersics")){
		i3_sersics_parameter_set * P = (i3_sersics_parameter_set*) p;
		if (!(P->x0==P->x0)) {return true;}
		if (!(P->y0==P->y0)) {return true;}
		if (!(P->e1==P->e1)) {return true;}
		if (!(P->e2==P->e2)) {return true;}
		if (!(P->radius==P->radius)) {return true;}
		if (!(P->radius_ratio==P->radius_ratio)) {return true;}
		if (!(P->bulge_A==P->bulge_A)) {return true;}
		if (!(P->disc_A==P->disc_A)) {return true;}
		if (!(P->bulge_index==P->bulge_index)) {return true;}
		if (!(P->disc_index==P->disc_index)) {return true;}
		if (!(P->delta_e_bulge==P->delta_e_bulge)) {return true;}
		if (!(P->delta_theta_bulge==P->delta_theta_bulge)) {return true;}
	}
	if (!strcmp(model->name,"logsersics")){
		i3_logsersics_parameter_set * P = (i3_logsersics_parameter_set*) p;
		if (!(P->x0==P->x0)) {return true;}
		if (!(P->y0==P->y0)) {return true;}
		if (!(P->e1==P->e1)) {return true;}
		if (!(P->e2==P->e2)) {return true;}
		if (!(P->radius==P->radius)) {return true;}
		if (!(P->radius_ratio==P->radius_ratio)) {return true;}
		if (!(P->bulge_A==P->bulge_A)) {return true;}
		if (!(P->disc_A==P->disc_A)) {return true;}
		if (!(P->bulge_index==P->bulge_index)) {return true;}
		if (!(P->disc_index==P->disc_index)) {return true;}
		if (!(P->delta_e_bulge==P->delta_e_bulge)) {return true;}
		if (!(P->delta_theta_bulge==P->delta_theta_bulge)) {return true;}
	}
	if (!strcmp(model->name,"quadratic_test")){
		i3_quadratic_test_parameter_set * P = (i3_quadratic_test_parameter_set*) p;
		if (!(P->x1==P->x1)) {return true;}
		if (!(P->x2==P->x2)) {return true;}
		if (!(P->x3==P->x3)) {return true;}
		if (!(P->x4==P->x4)) {return true;}
	}
	if (!strcmp(model->name,"multisersics")){
		i3_multisersics_parameter_set * P = (i3_multisersics_parameter_set*) p;
		if (!(P->ra_as==P->ra_as)) {return true;}
		if (!(P->dec_as==P->dec_as)) {return true;}
		if (!(P->e1==P->e1)) {return true;}
		if (!(P->e2==P->e2)) {return true;}
		if (!(P->radius==P->radius)) {return true;}
		if (!(P->radius_ratio==P->radius_ratio)) {return true;}
		if (!(P->bulge_A==P->bulge_A)) {return true;}
		if (!(P->disc_A==P->disc_A)) {return true;}
		if (!(P->bulge_index==P->bulge_index)) {return true;}
		if (!(P->disc_index==P->disc_index)) {return true;}
	}
	if (!strcmp(model->name,"multixray")){
		i3_multixray_parameter_set * P = (i3_multixray_parameter_set*) p;
		if (!(P->alpha_arcsec==P->alpha_arcsec)) {return true;}
		if (!(P->delta_arcsec==P->delta_arcsec)) {return true;}
		if (!(P->e1==P->e1)) {return true;}
		if (!(P->e2==P->e2)) {return true;}
		if (!(P->rc==P->rc)) {return true;}
		if (!(P->A==P->A)) {return true;}
		if (!(P->beta==P->beta)) {return true;}
	}
	return false;
}

void i3_model_option_fixes(i3_model * model, i3_options * options){
	if (!strcmp(model->name,"sersics")){
		model->param_fixed[0]=options->sersics_x0_fixed;
		model->param_fixed[1]=options->sersics_y0_fixed;
		model->param_fixed[2]=options->sersics_e1_fixed;
		model->param_fixed[3]=options->sersics_e2_fixed;
		model->param_fixed[4]=options->sersics_radius_fixed;
		model->param_fixed[5]=options->sersics_radius_ratio_fixed;
		model->param_fixed[6]=options->sersics_bulge_A_fixed;
		model->param_fixed[7]=options->sersics_disc_A_fixed;
		model->param_fixed[8]=options->sersics_bulge_index_fixed;
		model->param_fixed[9]=options->sersics_disc_index_fixed;
		model->param_fixed[10]=options->sersics_delta_e_bulge_fixed;
		model->param_fixed[11]=options->sersics_delta_theta_bulge_fixed;
}
	else if (!strcmp(model->name,"logsersics")){
		model->param_fixed[0]=options->logsersics_x0_fixed;
		model->param_fixed[1]=options->logsersics_y0_fixed;
		model->param_fixed[2]=options->logsersics_e1_fixed;
		model->param_fixed[3]=options->logsersics_e2_fixed;
		model->param_fixed[4]=options->logsersics_radius_fixed;
		model->param_fixed[5]=options->logsersics_radius_ratio_fixed;
		model->param_fixed[6]=options->logsersics_bulge_A_fixed;
		model->param_fixed[7]=options->logsersics_disc_A_fixed;
		model->param_fixed[8]=options->logsersics_bulge_index_fixed;
		model->param_fixed[9]=options->logsersics_disc_index_fixed;
		model->param_fixed[10]=options->logsersics_delta_e_bulge_fixed;
		model->param_fixed[11]=options->logsersics_delta_theta_bulge_fixed;
}
	else if (!strcmp(model->name,"quadratic_test")){
		model->param_fixed[0]=options->quadratic_test_x1_fixed;
		model->param_fixed[1]=options->quadratic_test_x2_fixed;
		model->param_fixed[2]=options->quadratic_test_x3_fixed;
		model->param_fixed[3]=options->quadratic_test_x4_fixed;
		model->param_fixed[4]=options->quadratic_test_x5_fixed;
}
	else if (!strcmp(model->name,"multisersics")){
		model->param_fixed[0]=options->multisersics_ra_as_fixed;
		model->param_fixed[1]=options->multisersics_dec_as_fixed;
		model->param_fixed[2]=options->multisersics_e1_fixed;
		model->param_fixed[3]=options->multisersics_e2_fixed;
		model->param_fixed[4]=options->multisersics_radius_fixed;
		model->param_fixed[5]=options->multisersics_radius_ratio_fixed;
		model->param_fixed[6]=options->multisersics_bulge_A_fixed;
		model->param_fixed[7]=options->multisersics_disc_A_fixed;
		model->param_fixed[8]=options->multisersics_bulge_index_fixed;
		model->param_fixed[9]=options->multisersics_disc_index_fixed;
}
	else if (!strcmp(model->name,"multixray")){
		model->param_fixed[0]=options->multixray_alpha_arcsec_fixed;
		model->param_fixed[1]=options->multixray_delta_arcsec_fixed;
		model->param_fixed[2]=options->multixray_e1_fixed;
		model->param_fixed[3]=options->multixray_e2_fixed;
		model->param_fixed[4]=options->multixray_rc_fixed;
		model->param_fixed[5]=options->multixray_A_fixed;
		model->param_fixed[6]=options->multixray_beta_fixed;
		model->param_fixed[7]=options->multixray_background_fixed;
}
}

int i3_model_set_parameter_by_name(i3_model * model, i3_parameter_set * p, char * name, double value)
{
		if (!strcmp(model->name,"sersics")){
		i3_sersics_parameter_set * P = (i3_sersics_parameter_set*) p;
		if (0==strcmp(name,"x0")) {P->x0=(float)value; return 0;}
		if (0==strcmp(name,"y0")) {P->y0=(float)value; return 0;}
		if (0==strcmp(name,"e1")) {P->e1=(float)value; return 0;}
		if (0==strcmp(name,"e2")) {P->e2=(float)value; return 0;}
		if (0==strcmp(name,"radius")) {P->radius=(float)value; return 0;}
		if (0==strcmp(name,"radius_ratio")) {P->radius_ratio=(float)value; return 0;}
		if (0==strcmp(name,"bulge_A")) {P->bulge_A=(float)value; return 0;}
		if (0==strcmp(name,"disc_A")) {P->disc_A=(float)value; return 0;}
		if (0==strcmp(name,"bulge_index")) {P->bulge_index=(float)value; return 0;}
		if (0==strcmp(name,"disc_index")) {P->disc_index=(float)value; return 0;}
		if (0==strcmp(name,"delta_e_bulge")) {P->delta_e_bulge=(float)value; return 0;}
		if (0==strcmp(name,"delta_theta_bulge")) {P->delta_theta_bulge=(float)value; return 0;}
}
	if (!strcmp(model->name,"logsersics")){
		i3_logsersics_parameter_set * P = (i3_logsersics_parameter_set*) p;
		if (0==strcmp(name,"x0")) {P->x0=(float)value; return 0;}
		if (0==strcmp(name,"y0")) {P->y0=(float)value; return 0;}
		if (0==strcmp(name,"e1")) {P->e1=(float)value; return 0;}
		if (0==strcmp(name,"e2")) {P->e2=(float)value; return 0;}
		if (0==strcmp(name,"radius")) {P->radius=(float)value; return 0;}
		if (0==strcmp(name,"radius_ratio")) {P->radius_ratio=(float)value; return 0;}
		if (0==strcmp(name,"bulge_A")) {P->bulge_A=(float)value; return 0;}
		if (0==strcmp(name,"disc_A")) {P->disc_A=(float)value; return 0;}
		if (0==strcmp(name,"bulge_index")) {P->bulge_index=(float)value; return 0;}
		if (0==strcmp(name,"disc_index")) {P->disc_index=(float)value; return 0;}
		if (0==strcmp(name,"delta_e_bulge")) {P->delta_e_bulge=(float)value; return 0;}
		if (0==strcmp(name,"delta_theta_bulge")) {P->delta_theta_bulge=(float)value; return 0;}
}
	if (!strcmp(model->name,"quadratic_test")){
		i3_quadratic_test_parameter_set * P = (i3_quadratic_test_parameter_set*) p;
		if (0==strcmp(name,"x1")) {P->x1=(float)value; return 0;}
		if (0==strcmp(name,"x2")) {P->x2=(float)value; return 0;}
		if (0==strcmp(name,"x3")) {P->x3=(float)value; return 0;}
		if (0==strcmp(name,"x4")) {P->x4=(float)value; return 0;}
		if (0==strcmp(name,"x5")) {fprintf(stderr, "Cannot set array params by name yet.\n"); return 1;}
}
	if (!strcmp(model->name,"multisersics")){
		i3_multisersics_parameter_set * P = (i3_multisersics_parameter_set*) p;
		if (0==strcmp(name,"ra_as")) {P->ra_as=(float)value; return 0;}
		if (0==strcmp(name,"dec_as")) {P->dec_as=(float)value; return 0;}
		if (0==strcmp(name,"e1")) {P->e1=(float)value; return 0;}
		if (0==strcmp(name,"e2")) {P->e2=(float)value; return 0;}
		if (0==strcmp(name,"radius")) {P->radius=(float)value; return 0;}
		if (0==strcmp(name,"radius_ratio")) {P->radius_ratio=(float)value; return 0;}
		if (0==strcmp(name,"bulge_A")) {P->bulge_A=(float)value; return 0;}
		if (0==strcmp(name,"disc_A")) {P->disc_A=(float)value; return 0;}
		if (0==strcmp(name,"bulge_index")) {P->bulge_index=(float)value; return 0;}
		if (0==strcmp(name,"disc_index")) {P->disc_index=(float)value; return 0;}
}
	if (!strcmp(model->name,"multixray")){
		i3_multixray_parameter_set * P = (i3_multixray_parameter_set*) p;
		if (0==strcmp(name,"alpha_arcsec")) {P->alpha_arcsec=(float)value; return 0;}
		if (0==strcmp(name,"delta_arcsec")) {P->delta_arcsec=(float)value; return 0;}
		if (0==strcmp(name,"e1")) {P->e1=(float)value; return 0;}
		if (0==strcmp(name,"e2")) {P->e2=(float)value; return 0;}
		if (0==strcmp(name,"rc")) {P->rc=(float)value; return 0;}
		if (0==strcmp(name,"A")) {P->A=(float)value; return 0;}
		if (0==strcmp(name,"beta")) {P->beta=(float)value; return 0;}
		if (0==strcmp(name,"background")) {fprintf(stderr, "Cannot set array params by name yet.\n"); return 1;}
}
	return 1;
}

i3_parameter_set * i3_model_option_starts(i3_model * model, i3_options * options){
	i3_parameter_set * p = (i3_parameter_set*)malloc(model->nbytes);
	if (!strcmp(model->name,"sersics")){
		i3_sersics_parameter_set * P = (i3_sersics_parameter_set*) p;
		P->x0=options->sersics_x0_start;
		P->y0=options->sersics_y0_start;
		P->e1=options->sersics_e1_start;
		P->e2=options->sersics_e2_start;
		P->radius=options->sersics_radius_start;
		P->radius_ratio=options->sersics_radius_ratio_start;
		P->bulge_A=options->sersics_bulge_A_start;
		P->disc_A=options->sersics_disc_A_start;
		P->bulge_index=options->sersics_bulge_index_start;
		P->disc_index=options->sersics_disc_index_start;
		P->delta_e_bulge=options->sersics_delta_e_bulge_start;
		P->delta_theta_bulge=options->sersics_delta_theta_bulge_start;
}
	else if (!strcmp(model->name,"logsersics")){
		i3_logsersics_parameter_set * P = (i3_logsersics_parameter_set*) p;
		P->x0=options->logsersics_x0_start;
		P->y0=options->logsersics_y0_start;
		P->e1=options->logsersics_e1_start;
		P->e2=options->logsersics_e2_start;
		P->radius=options->logsersics_radius_start;
		P->radius_ratio=options->logsersics_radius_ratio_start;
		P->bulge_A=options->logsersics_bulge_A_start;
		P->disc_A=options->logsersics_disc_A_start;
		P->bulge_index=options->logsersics_bulge_index_start;
		P->disc_index=options->logsersics_disc_index_start;
		P->delta_e_bulge=options->logsersics_delta_e_bulge_start;
		P->delta_theta_bulge=options->logsersics_delta_theta_bulge_start;
}
	else if (!strcmp(model->name,"quadratic_test")){
		i3_quadratic_test_parameter_set * P = (i3_quadratic_test_parameter_set*) p;
		P->x1=options->quadratic_test_x1_start;
		P->x2=options->quadratic_test_x2_start;
		P->x3=options->quadratic_test_x3_start;
		P->x4=options->quadratic_test_x4_start;
		for (int ip=0; ip<3; ip++) P->x5[ip]=options->quadratic_test_x5_start[ip];
}
	else if (!strcmp(model->name,"multisersics")){
		i3_multisersics_parameter_set * P = (i3_multisersics_parameter_set*) p;
		P->ra_as=options->multisersics_ra_as_start;
		P->dec_as=options->multisersics_dec_as_start;
		P->e1=options->multisersics_e1_start;
		P->e2=options->multisersics_e2_start;
		P->radius=options->multisersics_radius_start;
		P->radius_ratio=options->multisersics_radius_ratio_start;
		P->bulge_A=options->multisersics_bulge_A_start;
		P->disc_A=options->multisersics_disc_A_start;
		P->bulge_index=options->multisersics_bulge_index_start;
		P->disc_index=options->multisersics_disc_index_start;
}
	else if (!strcmp(model->name,"multixray")){
		i3_multixray_parameter_set * P = (i3_multixray_parameter_set*) p;
		P->alpha_arcsec=options->multixray_alpha_arcsec_start;
		P->delta_arcsec=options->multixray_delta_arcsec_start;
		P->e1=options->multixray_e1_start;
		P->e2=options->multixray_e2_start;
		P->rc=options->multixray_rc_start;
		P->A=options->multixray_A_start;
		P->beta=options->multixray_beta_start;
		for (int ip=0; ip<30; ip++) P->background[ip]=options->multixray_background_start[ip];
}
	return p;
}

i3_parameter_set * i3_model_option_min(i3_model * model, i3_options * options){
	i3_parameter_set * p = (i3_parameter_set*)malloc(model->nbytes);
	if (!strcmp(model->name,"sersics")){
		i3_sersics_parameter_set * P = (i3_sersics_parameter_set*) p;
		P->x0=options->sersics_x0_min;
		P->y0=options->sersics_y0_min;
		P->e1=options->sersics_e1_min;
		P->e2=options->sersics_e2_min;
		P->radius=options->sersics_radius_min;
		P->radius_ratio=options->sersics_radius_ratio_min;
		P->bulge_A=options->sersics_bulge_A_min;
		P->disc_A=options->sersics_disc_A_min;
		P->bulge_index=options->sersics_bulge_index_min;
		P->disc_index=options->sersics_disc_index_min;
		P->delta_e_bulge=options->sersics_delta_e_bulge_min;
		P->delta_theta_bulge=options->sersics_delta_theta_bulge_min;
}
	else if (!strcmp(model->name,"logsersics")){
		i3_logsersics_parameter_set * P = (i3_logsersics_parameter_set*) p;
		P->x0=options->logsersics_x0_min;
		P->y0=options->logsersics_y0_min;
		P->e1=options->logsersics_e1_min;
		P->e2=options->logsersics_e2_min;
		P->radius=options->logsersics_radius_min;
		P->radius_ratio=options->logsersics_radius_ratio_min;
		P->bulge_A=options->logsersics_bulge_A_min;
		P->disc_A=options->logsersics_disc_A_min;
		P->bulge_index=options->logsersics_bulge_index_min;
		P->disc_index=options->logsersics_disc_index_min;
		P->delta_e_bulge=options->logsersics_delta_e_bulge_min;
		P->delta_theta_bulge=options->logsersics_delta_theta_bulge_min;
}
	else if (!strcmp(model->name,"quadratic_test")){
		i3_quadratic_test_parameter_set * P = (i3_quadratic_test_parameter_set*) p;
		P->x1=options->quadratic_test_x1_min;
		P->x2=options->quadratic_test_x2_min;
		P->x3=options->quadratic_test_x3_min;
		P->x4=options->quadratic_test_x4_min;
		for (int ip=0; ip<3; ip++) P->x5[ip]=options->quadratic_test_x5_min[ip];
}
	else if (!strcmp(model->name,"multisersics")){
		i3_multisersics_parameter_set * P = (i3_multisersics_parameter_set*) p;
		P->ra_as=options->multisersics_ra_as_min;
		P->dec_as=options->multisersics_dec_as_min;
		P->e1=options->multisersics_e1_min;
		P->e2=options->multisersics_e2_min;
		P->radius=options->multisersics_radius_min;
		P->radius_ratio=options->multisersics_radius_ratio_min;
		P->bulge_A=options->multisersics_bulge_A_min;
		P->disc_A=options->multisersics_disc_A_min;
		P->bulge_index=options->multisersics_bulge_index_min;
		P->disc_index=options->multisersics_disc_index_min;
}
	else if (!strcmp(model->name,"multixray")){
		i3_multixray_parameter_set * P = (i3_multixray_parameter_set*) p;
		P->alpha_arcsec=options->multixray_alpha_arcsec_min;
		P->delta_arcsec=options->multixray_delta_arcsec_min;
		P->e1=options->multixray_e1_min;
		P->e2=options->multixray_e2_min;
		P->rc=options->multixray_rc_min;
		P->A=options->multixray_A_min;
		P->beta=options->multixray_beta_min;
		for (int ip=0; ip<30; ip++) P->background[ip]=options->multixray_background_min[ip];
}
	return p;
}

i3_parameter_set * i3_model_option_max(i3_model * model, i3_options * options){
	i3_parameter_set * p = (i3_parameter_set*)malloc(model->nbytes);
	if (!strcmp(model->name,"sersics")){
		i3_sersics_parameter_set * P = (i3_sersics_parameter_set*) p;
		P->x0=options->sersics_x0_max;
		P->y0=options->sersics_y0_max;
		P->e1=options->sersics_e1_max;
		P->e2=options->sersics_e2_max;
		P->radius=options->sersics_radius_max;
		P->radius_ratio=options->sersics_radius_ratio_max;
		P->bulge_A=options->sersics_bulge_A_max;
		P->disc_A=options->sersics_disc_A_max;
		P->bulge_index=options->sersics_bulge_index_max;
		P->disc_index=options->sersics_disc_index_max;
		P->delta_e_bulge=options->sersics_delta_e_bulge_max;
		P->delta_theta_bulge=options->sersics_delta_theta_bulge_max;
}
	else if (!strcmp(model->name,"logsersics")){
		i3_logsersics_parameter_set * P = (i3_logsersics_parameter_set*) p;
		P->x0=options->logsersics_x0_max;
		P->y0=options->logsersics_y0_max;
		P->e1=options->logsersics_e1_max;
		P->e2=options->logsersics_e2_max;
		P->radius=options->logsersics_radius_max;
		P->radius_ratio=options->logsersics_radius_ratio_max;
		P->bulge_A=options->logsersics_bulge_A_max;
		P->disc_A=options->logsersics_disc_A_max;
		P->bulge_index=options->logsersics_bulge_index_max;
		P->disc_index=options->logsersics_disc_index_max;
		P->delta_e_bulge=options->logsersics_delta_e_bulge_max;
		P->delta_theta_bulge=options->logsersics_delta_theta_bulge_max;
}
	else if (!strcmp(model->name,"quadratic_test")){
		i3_quadratic_test_parameter_set * P = (i3_quadratic_test_parameter_set*) p;
		P->x1=options->quadratic_test_x1_max;
		P->x2=options->quadratic_test_x2_max;
		P->x3=options->quadratic_test_x3_max;
		P->x4=options->quadratic_test_x4_max;
		for (int ip=0; ip<3; ip++) P->x5[ip]=options->quadratic_test_x5_max[ip];
}
	else if (!strcmp(model->name,"multisersics")){
		i3_multisersics_parameter_set * P = (i3_multisersics_parameter_set*) p;
		P->ra_as=options->multisersics_ra_as_max;
		P->dec_as=options->multisersics_dec_as_max;
		P->e1=options->multisersics_e1_max;
		P->e2=options->multisersics_e2_max;
		P->radius=options->multisersics_radius_max;
		P->radius_ratio=options->multisersics_radius_ratio_max;
		P->bulge_A=options->multisersics_bulge_A_max;
		P->disc_A=options->multisersics_disc_A_max;
		P->bulge_index=options->multisersics_bulge_index_max;
		P->disc_index=options->multisersics_disc_index_max;
}
	else if (!strcmp(model->name,"multixray")){
		i3_multixray_parameter_set * P = (i3_multixray_parameter_set*) p;
		P->alpha_arcsec=options->multixray_alpha_arcsec_max;
		P->delta_arcsec=options->multixray_delta_arcsec_max;
		P->e1=options->multixray_e1_max;
		P->e2=options->multixray_e2_max;
		P->rc=options->multixray_rc_max;
		P->A=options->multixray_A_max;
		P->beta=options->multixray_beta_max;
		for (int ip=0; ip<30; ip++) P->background[ip]=options->multixray_background_max[ip];
}
	return p;
}

i3_parameter_set * i3_model_option_range(i3_model * model, i3_options * options){
	i3_parameter_set * p = (i3_parameter_set*)malloc(model->nbytes);
	if (!strcmp(model->name,"sersics")){
		i3_sersics_parameter_set * P = (i3_sersics_parameter_set*) p;
		i3_sersics_parameter_set * range_min = (i3_sersics_parameter_set*) i3_model_option_min(model,options);
		i3_sersics_parameter_set * range_max = (i3_sersics_parameter_set*) i3_model_option_max(model,options);
		P->x0=range_max->x0-range_min->x0;
		P->y0=range_max->y0-range_min->y0;
		P->e1=range_max->e1-range_min->e1;
		P->e2=range_max->e2-range_min->e2;
		P->radius=range_max->radius-range_min->radius;
		P->radius_ratio=range_max->radius_ratio-range_min->radius_ratio;
		P->bulge_A=range_max->bulge_A-range_min->bulge_A;
		P->disc_A=range_max->disc_A-range_min->disc_A;
		P->bulge_index=range_max->bulge_index-range_min->bulge_index;
		P->disc_index=range_max->disc_index-range_min->disc_index;
		P->delta_e_bulge=range_max->delta_e_bulge-range_min->delta_e_bulge;
		P->delta_theta_bulge=range_max->delta_theta_bulge-range_min->delta_theta_bulge;
		free(range_min);
		free(range_max);
}
	else if (!strcmp(model->name,"logsersics")){
		i3_logsersics_parameter_set * P = (i3_logsersics_parameter_set*) p;
		i3_logsersics_parameter_set * range_min = (i3_logsersics_parameter_set*) i3_model_option_min(model,options);
		i3_logsersics_parameter_set * range_max = (i3_logsersics_parameter_set*) i3_model_option_max(model,options);
		P->x0=range_max->x0-range_min->x0;
		P->y0=range_max->y0-range_min->y0;
		P->e1=range_max->e1-range_min->e1;
		P->e2=range_max->e2-range_min->e2;
		P->radius=range_max->radius-range_min->radius;
		P->radius_ratio=range_max->radius_ratio-range_min->radius_ratio;
		P->bulge_A=range_max->bulge_A-range_min->bulge_A;
		P->disc_A=range_max->disc_A-range_min->disc_A;
		P->bulge_index=range_max->bulge_index-range_min->bulge_index;
		P->disc_index=range_max->disc_index-range_min->disc_index;
		P->delta_e_bulge=range_max->delta_e_bulge-range_min->delta_e_bulge;
		P->delta_theta_bulge=range_max->delta_theta_bulge-range_min->delta_theta_bulge;
		free(range_min);
		free(range_max);
}
	else if (!strcmp(model->name,"quadratic_test")){
		i3_quadratic_test_parameter_set * P = (i3_quadratic_test_parameter_set*) p;
		i3_quadratic_test_parameter_set * range_min = (i3_quadratic_test_parameter_set*) i3_model_option_min(model,options);
		i3_quadratic_test_parameter_set * range_max = (i3_quadratic_test_parameter_set*) i3_model_option_max(model,options);
		P->x1=range_max->x1-range_min->x1;
		P->x2=range_max->x2-range_min->x2;
		P->x3=range_max->x3-range_min->x3;
		P->x4=range_max->x4-range_min->x4;
		for (int ip=0; ip<3; ip++) P->x5[ip]=range_max->x5[ip]-range_min->x5[ip];
		free(range_min);
		free(range_max);
}
	else if (!strcmp(model->name,"multisersics")){
		i3_multisersics_parameter_set * P = (i3_multisersics_parameter_set*) p;
		i3_multisersics_parameter_set * range_min = (i3_multisersics_parameter_set*) i3_model_option_min(model,options);
		i3_multisersics_parameter_set * range_max = (i3_multisersics_parameter_set*) i3_model_option_max(model,options);
		P->ra_as=range_max->ra_as-range_min->ra_as;
		P->dec_as=range_max->dec_as-range_min->dec_as;
		P->e1=range_max->e1-range_min->e1;
		P->e2=range_max->e2-range_min->e2;
		P->radius=range_max->radius-range_min->radius;
		P->radius_ratio=range_max->radius_ratio-range_min->radius_ratio;
		P->bulge_A=range_max->bulge_A-range_min->bulge_A;
		P->disc_A=range_max->disc_A-range_min->disc_A;
		P->bulge_index=range_max->bulge_index-range_min->bulge_index;
		P->disc_index=range_max->disc_index-range_min->disc_index;
		free(range_min);
		free(range_max);
}
	else if (!strcmp(model->name,"multixray")){
		i3_multixray_parameter_set * P = (i3_multixray_parameter_set*) p;
		i3_multixray_parameter_set * range_min = (i3_multixray_parameter_set*) i3_model_option_min(model,options);
		i3_multixray_parameter_set * range_max = (i3_multixray_parameter_set*) i3_model_option_max(model,options);
		P->alpha_arcsec=range_max->alpha_arcsec-range_min->alpha_arcsec;
		P->delta_arcsec=range_max->delta_arcsec-range_min->delta_arcsec;
		P->e1=range_max->e1-range_min->e1;
		P->e2=range_max->e2-range_min->e2;
		P->rc=range_max->rc-range_min->rc;
		P->A=range_max->A-range_min->A;
		P->beta=range_max->beta-range_min->beta;
		for (int ip=0; ip<30; ip++) P->background[ip]=range_max->background[ip]-range_min->background[ip];
		free(range_min);
		free(range_max);
}
	return p;
}


void i3_model_posterior_derivative(i3_model * model, i3_image * model_image, i3_flt likelihood, i3_data_set * data_set, i3_parameter_set * parameters, i3_parameter_set * derivative){
	if (!model->posterior_derivative) I3_FATAL("Attempted to compute derivative of model, but not available",22);
	model->posterior_derivative(model_image,likelihood,data_set,parameters,derivative);
	
}




int i3_model_number_varied_params(i3_model * model){
	int n=0;
	for(int i=0;i<model->nparam;i++){
		if (!model->param_fixed[i]) n++;
	}
	return n;
}

int i3_model_number_varied_params_nonzero_width(i3_model * model, i3_parameter_set * width){
	int n=0;
	if (!strcmp(model->name,"sersics")){
		i3_sersics_parameter_set * W = (i3_sersics_parameter_set*) width;
		if (W->x0!=0 && !model->param_fixed[0]) n+=1;
		if (W->y0!=0 && !model->param_fixed[1]) n+=1;
		if (W->e1!=0 && !model->param_fixed[2]) n+=1;
		if (W->e2!=0 && !model->param_fixed[3]) n+=1;
		if (W->radius!=0 && !model->param_fixed[4]) n+=1;
		if (W->radius_ratio!=0 && !model->param_fixed[5]) n+=1;
		if (W->bulge_A!=0 && !model->param_fixed[6]) n+=1;
		if (W->disc_A!=0 && !model->param_fixed[7]) n+=1;
		if (W->bulge_index!=0 && !model->param_fixed[8]) n+=1;
		if (W->disc_index!=0 && !model->param_fixed[9]) n+=1;
		if (W->delta_e_bulge!=0 && !model->param_fixed[10]) n+=1;
		if (W->delta_theta_bulge!=0 && !model->param_fixed[11]) n+=1;
	}
	else if (!strcmp(model->name,"logsersics")){
		i3_logsersics_parameter_set * W = (i3_logsersics_parameter_set*) width;
		if (W->x0!=0 && !model->param_fixed[0]) n+=1;
		if (W->y0!=0 && !model->param_fixed[1]) n+=1;
		if (W->e1!=0 && !model->param_fixed[2]) n+=1;
		if (W->e2!=0 && !model->param_fixed[3]) n+=1;
		if (W->radius!=0 && !model->param_fixed[4]) n+=1;
		if (W->radius_ratio!=0 && !model->param_fixed[5]) n+=1;
		if (W->bulge_A!=0 && !model->param_fixed[6]) n+=1;
		if (W->disc_A!=0 && !model->param_fixed[7]) n+=1;
		if (W->bulge_index!=0 && !model->param_fixed[8]) n+=1;
		if (W->disc_index!=0 && !model->param_fixed[9]) n+=1;
		if (W->delta_e_bulge!=0 && !model->param_fixed[10]) n+=1;
		if (W->delta_theta_bulge!=0 && !model->param_fixed[11]) n+=1;
	}
	else if (!strcmp(model->name,"quadratic_test")){
		i3_quadratic_test_parameter_set * W = (i3_quadratic_test_parameter_set*) width;
		if (W->x1!=0 && !model->param_fixed[0]) n+=1;
		if (W->x2!=0 && !model->param_fixed[1]) n+=1;
		if (W->x3!=0 && !model->param_fixed[2]) n+=1;
		if (W->x4!=0 && !model->param_fixed[3]) n+=1;
		if (W->x5!=0 && !model->param_fixed[4]) n+=1;
	}
	else if (!strcmp(model->name,"multisersics")){
		i3_multisersics_parameter_set * W = (i3_multisersics_parameter_set*) width;
		if (W->ra_as!=0 && !model->param_fixed[0]) n+=1;
		if (W->dec_as!=0 && !model->param_fixed[1]) n+=1;
		if (W->e1!=0 && !model->param_fixed[2]) n+=1;
		if (W->e2!=0 && !model->param_fixed[3]) n+=1;
		if (W->radius!=0 && !model->param_fixed[4]) n+=1;
		if (W->radius_ratio!=0 && !model->param_fixed[5]) n+=1;
		if (W->bulge_A!=0 && !model->param_fixed[6]) n+=1;
		if (W->disc_A!=0 && !model->param_fixed[7]) n+=1;
		if (W->bulge_index!=0 && !model->param_fixed[8]) n+=1;
		if (W->disc_index!=0 && !model->param_fixed[9]) n+=1;
	}
	else if (!strcmp(model->name,"multixray")){
		i3_multixray_parameter_set * W = (i3_multixray_parameter_set*) width;
		if (W->alpha_arcsec!=0 && !model->param_fixed[0]) n+=1;
		if (W->delta_arcsec!=0 && !model->param_fixed[1]) n+=1;
		if (W->e1!=0 && !model->param_fixed[2]) n+=1;
		if (W->e2!=0 && !model->param_fixed[3]) n+=1;
		if (W->rc!=0 && !model->param_fixed[4]) n+=1;
		if (W->A!=0 && !model->param_fixed[5]) n+=1;
		if (W->beta!=0 && !model->param_fixed[6]) n+=1;
		if (W->background!=0 && !model->param_fixed[7]) n+=1;
	}
	return n;
}


void i3_model_copy_parameters(i3_model * model, i3_parameter_set * restrict dest, i3_parameter_set * restrict source){
	memcpy(dest,source,model->nbytes);
}

void i3_model_scale_parameters(i3_model * model, i3_parameter_set * p, i3_flt scale){
	if (!strcmp(model->name,"sersics")){
		i3_sersics_parameter_set * P = (i3_sersics_parameter_set*) p;
		P->x0*=scale;
		P->y0*=scale;
		P->e1*=scale;
		P->e2*=scale;
		P->radius*=scale;
		P->radius_ratio*=scale;
		P->bulge_A*=scale;
		P->disc_A*=scale;
		P->bulge_index*=scale;
		P->disc_index*=scale;
		P->delta_e_bulge*=scale;
		P->delta_theta_bulge*=scale;

	}
	else if (!strcmp(model->name,"logsersics")){
		i3_logsersics_parameter_set * P = (i3_logsersics_parameter_set*) p;
		P->x0*=scale;
		P->y0*=scale;
		P->e1*=scale;
		P->e2*=scale;
		P->radius*=scale;
		P->radius_ratio*=scale;
		P->bulge_A*=scale;
		P->disc_A*=scale;
		P->bulge_index*=scale;
		P->disc_index*=scale;
		P->delta_e_bulge*=scale;
		P->delta_theta_bulge*=scale;

	}
	else if (!strcmp(model->name,"quadratic_test")){
		i3_quadratic_test_parameter_set * P = (i3_quadratic_test_parameter_set*) p;
		P->x1*=scale;
		P->x2*=scale;
		P->x3*=scale;
		P->x4*=scale;
		for (int ip=0;ip<3;ip++) P->x5[ip]*=scale;

	}
	else if (!strcmp(model->name,"multisersics")){
		i3_multisersics_parameter_set * P = (i3_multisersics_parameter_set*) p;
		P->ra_as*=scale;
		P->dec_as*=scale;
		P->e1*=scale;
		P->e2*=scale;
		P->radius*=scale;
		P->radius_ratio*=scale;
		P->bulge_A*=scale;
		P->disc_A*=scale;
		P->bulge_index*=scale;
		P->disc_index*=scale;

	}
	else if (!strcmp(model->name,"multixray")){
		i3_multixray_parameter_set * P = (i3_multixray_parameter_set*) p;
		P->alpha_arcsec*=scale;
		P->delta_arcsec*=scale;
		P->e1*=scale;
		P->e2*=scale;
		P->rc*=scale;
		P->A*=scale;
		P->beta*=scale;
		for (int ip=0;ip<30;ip++) P->background[ip]*=scale;

	}
}

void i3_model_perturb(i3_model * model, i3_parameter_set * p, i3_flt scale)
{
		if (!strcmp(model->name,"sersics")){
		i3_sersics_parameter_set * P = (i3_sersics_parameter_set*) p;
		i3_sersics_parameter_set * min_p = (i3_sersics_parameter_set*)model->min;
		i3_sersics_parameter_set * max_p = (i3_sersics_parameter_set*)model->max;
		i3_flt epsilon = 1e-6;
		 if (!model->param_fixed[0]) {
			i3_flt r = (max_p->x0 - min_p->x0);
			P->x0 += scale * r * i3_random_normal();
			if (P->x0<min_p->x0) P->x0 = min_p->x0+r*epsilon;
			if (P->x0>max_p->x0) P->x0 = max_p->x0-r*epsilon;
		}
		 if (!model->param_fixed[1]) {
			i3_flt r = (max_p->y0 - min_p->y0);
			P->y0 += scale * r * i3_random_normal();
			if (P->y0<min_p->y0) P->y0 = min_p->y0+r*epsilon;
			if (P->y0>max_p->y0) P->y0 = max_p->y0-r*epsilon;
		}
		 if (!model->param_fixed[2]) {
			i3_flt r = (max_p->e1 - min_p->e1);
			P->e1 += scale * r * i3_random_normal();
			if (P->e1<min_p->e1) P->e1 = min_p->e1+r*epsilon;
			if (P->e1>max_p->e1) P->e1 = max_p->e1-r*epsilon;
		}
		 if (!model->param_fixed[3]) {
			i3_flt r = (max_p->e2 - min_p->e2);
			P->e2 += scale * r * i3_random_normal();
			if (P->e2<min_p->e2) P->e2 = min_p->e2+r*epsilon;
			if (P->e2>max_p->e2) P->e2 = max_p->e2-r*epsilon;
		}
		 if (!model->param_fixed[4]) {
			i3_flt r = (max_p->radius - min_p->radius);
			P->radius += scale * r * i3_random_normal();
			if (P->radius<min_p->radius) P->radius = min_p->radius+r*epsilon;
			if (P->radius>max_p->radius) P->radius = max_p->radius-r*epsilon;
		}
		 if (!model->param_fixed[5]) {
			i3_flt r = (max_p->radius_ratio - min_p->radius_ratio);
			P->radius_ratio += scale * r * i3_random_normal();
			if (P->radius_ratio<min_p->radius_ratio) P->radius_ratio = min_p->radius_ratio+r*epsilon;
			if (P->radius_ratio>max_p->radius_ratio) P->radius_ratio = max_p->radius_ratio-r*epsilon;
		}
		 if (!model->param_fixed[6]) {
			i3_flt r = (max_p->bulge_A - min_p->bulge_A);
			P->bulge_A += scale * r * i3_random_normal();
			if (P->bulge_A<min_p->bulge_A) P->bulge_A = min_p->bulge_A+r*epsilon;
			if (P->bulge_A>max_p->bulge_A) P->bulge_A = max_p->bulge_A-r*epsilon;
		}
		 if (!model->param_fixed[7]) {
			i3_flt r = (max_p->disc_A - min_p->disc_A);
			P->disc_A += scale * r * i3_random_normal();
			if (P->disc_A<min_p->disc_A) P->disc_A = min_p->disc_A+r*epsilon;
			if (P->disc_A>max_p->disc_A) P->disc_A = max_p->disc_A-r*epsilon;
		}
		 if (!model->param_fixed[8]) {
			i3_flt r = (max_p->bulge_index - min_p->bulge_index);
			P->bulge_index += scale * r * i3_random_normal();
			if (P->bulge_index<min_p->bulge_index) P->bulge_index = min_p->bulge_index+r*epsilon;
			if (P->bulge_index>max_p->bulge_index) P->bulge_index = max_p->bulge_index-r*epsilon;
		}
		 if (!model->param_fixed[9]) {
			i3_flt r = (max_p->disc_index - min_p->disc_index);
			P->disc_index += scale * r * i3_random_normal();
			if (P->disc_index<min_p->disc_index) P->disc_index = min_p->disc_index+r*epsilon;
			if (P->disc_index>max_p->disc_index) P->disc_index = max_p->disc_index-r*epsilon;
		}
		 if (!model->param_fixed[10]) {
			i3_flt r = (max_p->delta_e_bulge - min_p->delta_e_bulge);
			P->delta_e_bulge += scale * r * i3_random_normal();
			if (P->delta_e_bulge<min_p->delta_e_bulge) P->delta_e_bulge = min_p->delta_e_bulge+r*epsilon;
			if (P->delta_e_bulge>max_p->delta_e_bulge) P->delta_e_bulge = max_p->delta_e_bulge-r*epsilon;
		}
		 if (!model->param_fixed[11]) {
			i3_flt r = (max_p->delta_theta_bulge - min_p->delta_theta_bulge);
			P->delta_theta_bulge += scale * r * i3_random_normal();
			if (P->delta_theta_bulge<min_p->delta_theta_bulge) P->delta_theta_bulge = min_p->delta_theta_bulge+r*epsilon;
			if (P->delta_theta_bulge>max_p->delta_theta_bulge) P->delta_theta_bulge = max_p->delta_theta_bulge-r*epsilon;
		}
	}
	if (!strcmp(model->name,"logsersics")){
		i3_logsersics_parameter_set * P = (i3_logsersics_parameter_set*) p;
		i3_logsersics_parameter_set * min_p = (i3_logsersics_parameter_set*)model->min;
		i3_logsersics_parameter_set * max_p = (i3_logsersics_parameter_set*)model->max;
		i3_flt epsilon = 1e-6;
		 if (!model->param_fixed[0]) {
			i3_flt r = (max_p->x0 - min_p->x0);
			P->x0 += scale * r * i3_random_normal();
			if (P->x0<min_p->x0) P->x0 = min_p->x0+r*epsilon;
			if (P->x0>max_p->x0) P->x0 = max_p->x0-r*epsilon;
		}
		 if (!model->param_fixed[1]) {
			i3_flt r = (max_p->y0 - min_p->y0);
			P->y0 += scale * r * i3_random_normal();
			if (P->y0<min_p->y0) P->y0 = min_p->y0+r*epsilon;
			if (P->y0>max_p->y0) P->y0 = max_p->y0-r*epsilon;
		}
		 if (!model->param_fixed[2]) {
			i3_flt r = (max_p->e1 - min_p->e1);
			P->e1 += scale * r * i3_random_normal();
			if (P->e1<min_p->e1) P->e1 = min_p->e1+r*epsilon;
			if (P->e1>max_p->e1) P->e1 = max_p->e1-r*epsilon;
		}
		 if (!model->param_fixed[3]) {
			i3_flt r = (max_p->e2 - min_p->e2);
			P->e2 += scale * r * i3_random_normal();
			if (P->e2<min_p->e2) P->e2 = min_p->e2+r*epsilon;
			if (P->e2>max_p->e2) P->e2 = max_p->e2-r*epsilon;
		}
		 if (!model->param_fixed[4]) {
			i3_flt r = (max_p->radius - min_p->radius);
			P->radius += scale * r * i3_random_normal();
			if (P->radius<min_p->radius) P->radius = min_p->radius+r*epsilon;
			if (P->radius>max_p->radius) P->radius = max_p->radius-r*epsilon;
		}
		 if (!model->param_fixed[5]) {
			i3_flt r = (max_p->radius_ratio - min_p->radius_ratio);
			P->radius_ratio += scale * r * i3_random_normal();
			if (P->radius_ratio<min_p->radius_ratio) P->radius_ratio = min_p->radius_ratio+r*epsilon;
			if (P->radius_ratio>max_p->radius_ratio) P->radius_ratio = max_p->radius_ratio-r*epsilon;
		}
		 if (!model->param_fixed[6]) {
			i3_flt r = (max_p->bulge_A - min_p->bulge_A);
			P->bulge_A += scale * r * i3_random_normal();
			if (P->bulge_A<min_p->bulge_A) P->bulge_A = min_p->bulge_A+r*epsilon;
			if (P->bulge_A>max_p->bulge_A) P->bulge_A = max_p->bulge_A-r*epsilon;
		}
		 if (!model->param_fixed[7]) {
			i3_flt r = (max_p->disc_A - min_p->disc_A);
			P->disc_A += scale * r * i3_random_normal();
			if (P->disc_A<min_p->disc_A) P->disc_A = min_p->disc_A+r*epsilon;
			if (P->disc_A>max_p->disc_A) P->disc_A = max_p->disc_A-r*epsilon;
		}
		 if (!model->param_fixed[8]) {
			i3_flt r = (max_p->bulge_index - min_p->bulge_index);
			P->bulge_index += scale * r * i3_random_normal();
			if (P->bulge_index<min_p->bulge_index) P->bulge_index = min_p->bulge_index+r*epsilon;
			if (P->bulge_index>max_p->bulge_index) P->bulge_index = max_p->bulge_index-r*epsilon;
		}
		 if (!model->param_fixed[9]) {
			i3_flt r = (max_p->disc_index - min_p->disc_index);
			P->disc_index += scale * r * i3_random_normal();
			if (P->disc_index<min_p->disc_index) P->disc_index = min_p->disc_index+r*epsilon;
			if (P->disc_index>max_p->disc_index) P->disc_index = max_p->disc_index-r*epsilon;
		}
		 if (!model->param_fixed[10]) {
			i3_flt r = (max_p->delta_e_bulge - min_p->delta_e_bulge);
			P->delta_e_bulge += scale * r * i3_random_normal();
			if (P->delta_e_bulge<min_p->delta_e_bulge) P->delta_e_bulge = min_p->delta_e_bulge+r*epsilon;
			if (P->delta_e_bulge>max_p->delta_e_bulge) P->delta_e_bulge = max_p->delta_e_bulge-r*epsilon;
		}
		 if (!model->param_fixed[11]) {
			i3_flt r = (max_p->delta_theta_bulge - min_p->delta_theta_bulge);
			P->delta_theta_bulge += scale * r * i3_random_normal();
			if (P->delta_theta_bulge<min_p->delta_theta_bulge) P->delta_theta_bulge = min_p->delta_theta_bulge+r*epsilon;
			if (P->delta_theta_bulge>max_p->delta_theta_bulge) P->delta_theta_bulge = max_p->delta_theta_bulge-r*epsilon;
		}
	}
	if (!strcmp(model->name,"quadratic_test")){
		i3_quadratic_test_parameter_set * P = (i3_quadratic_test_parameter_set*) p;
		i3_quadratic_test_parameter_set * min_p = (i3_quadratic_test_parameter_set*)model->min;
		i3_quadratic_test_parameter_set * max_p = (i3_quadratic_test_parameter_set*)model->max;
		i3_flt epsilon = 1e-6;
		 if (!model->param_fixed[0]) {
			i3_flt r = (max_p->x1 - min_p->x1);
			P->x1 += scale * r * i3_random_normal();
			if (P->x1<min_p->x1) P->x1 = min_p->x1+r*epsilon;
			if (P->x1>max_p->x1) P->x1 = max_p->x1-r*epsilon;
		}
		 if (!model->param_fixed[1]) {
			i3_flt r = (max_p->x2 - min_p->x2);
			P->x2 += scale * r * i3_random_normal();
			if (P->x2<min_p->x2) P->x2 = min_p->x2+r*epsilon;
			if (P->x2>max_p->x2) P->x2 = max_p->x2-r*epsilon;
		}
		 if (!model->param_fixed[2]) {
			i3_flt r = (max_p->x3 - min_p->x3);
			P->x3 += scale * r * i3_random_normal();
			if (P->x3<min_p->x3) P->x3 = min_p->x3+r*epsilon;
			if (P->x3>max_p->x3) P->x3 = max_p->x3-r*epsilon;
		}
		 if (!model->param_fixed[3]) {
			i3_flt r = (max_p->x4 - min_p->x4);
			P->x4 += scale * r * i3_random_normal();
			if (P->x4<min_p->x4) P->x4 = min_p->x4+r*epsilon;
			if (P->x4>max_p->x4) P->x4 = max_p->x4-r*epsilon;
		}
		 if (!model->param_fixed[4]) {
			 for (int ip=0; ip<3; ip++){
				i3_flt r = (max_p->x5[ip] - min_p->x5[ip]);
				P->x5[ip] += scale * r * i3_random_normal();
				if (P->x5[ip]<min_p->x5[ip]) P->x5[ip] = min_p->x5[ip]+r*epsilon;
				if (P->x5[ip]>max_p->x5[ip]) P->x5[ip] = max_p->x5[ip]-r*epsilon;
			}
		}
	}
	if (!strcmp(model->name,"multisersics")){
		i3_multisersics_parameter_set * P = (i3_multisersics_parameter_set*) p;
		i3_multisersics_parameter_set * min_p = (i3_multisersics_parameter_set*)model->min;
		i3_multisersics_parameter_set * max_p = (i3_multisersics_parameter_set*)model->max;
		i3_flt epsilon = 1e-6;
		 if (!model->param_fixed[0]) {
			i3_flt r = (max_p->ra_as - min_p->ra_as);
			P->ra_as += scale * r * i3_random_normal();
			if (P->ra_as<min_p->ra_as) P->ra_as = min_p->ra_as+r*epsilon;
			if (P->ra_as>max_p->ra_as) P->ra_as = max_p->ra_as-r*epsilon;
		}
		 if (!model->param_fixed[1]) {
			i3_flt r = (max_p->dec_as - min_p->dec_as);
			P->dec_as += scale * r * i3_random_normal();
			if (P->dec_as<min_p->dec_as) P->dec_as = min_p->dec_as+r*epsilon;
			if (P->dec_as>max_p->dec_as) P->dec_as = max_p->dec_as-r*epsilon;
		}
		 if (!model->param_fixed[2]) {
			i3_flt r = (max_p->e1 - min_p->e1);
			P->e1 += scale * r * i3_random_normal();
			if (P->e1<min_p->e1) P->e1 = min_p->e1+r*epsilon;
			if (P->e1>max_p->e1) P->e1 = max_p->e1-r*epsilon;
		}
		 if (!model->param_fixed[3]) {
			i3_flt r = (max_p->e2 - min_p->e2);
			P->e2 += scale * r * i3_random_normal();
			if (P->e2<min_p->e2) P->e2 = min_p->e2+r*epsilon;
			if (P->e2>max_p->e2) P->e2 = max_p->e2-r*epsilon;
		}
		 if (!model->param_fixed[4]) {
			i3_flt r = (max_p->radius - min_p->radius);
			P->radius += scale * r * i3_random_normal();
			if (P->radius<min_p->radius) P->radius = min_p->radius+r*epsilon;
			if (P->radius>max_p->radius) P->radius = max_p->radius-r*epsilon;
		}
		 if (!model->param_fixed[5]) {
			i3_flt r = (max_p->radius_ratio - min_p->radius_ratio);
			P->radius_ratio += scale * r * i3_random_normal();
			if (P->radius_ratio<min_p->radius_ratio) P->radius_ratio = min_p->radius_ratio+r*epsilon;
			if (P->radius_ratio>max_p->radius_ratio) P->radius_ratio = max_p->radius_ratio-r*epsilon;
		}
		 if (!model->param_fixed[6]) {
			i3_flt r = (max_p->bulge_A - min_p->bulge_A);
			P->bulge_A += scale * r * i3_random_normal();
			if (P->bulge_A<min_p->bulge_A) P->bulge_A = min_p->bulge_A+r*epsilon;
			if (P->bulge_A>max_p->bulge_A) P->bulge_A = max_p->bulge_A-r*epsilon;
		}
		 if (!model->param_fixed[7]) {
			i3_flt r = (max_p->disc_A - min_p->disc_A);
			P->disc_A += scale * r * i3_random_normal();
			if (P->disc_A<min_p->disc_A) P->disc_A = min_p->disc_A+r*epsilon;
			if (P->disc_A>max_p->disc_A) P->disc_A = max_p->disc_A-r*epsilon;
		}
		 if (!model->param_fixed[8]) {
			i3_flt r = (max_p->bulge_index - min_p->bulge_index);
			P->bulge_index += scale * r * i3_random_normal();
			if (P->bulge_index<min_p->bulge_index) P->bulge_index = min_p->bulge_index+r*epsilon;
			if (P->bulge_index>max_p->bulge_index) P->bulge_index = max_p->bulge_index-r*epsilon;
		}
		 if (!model->param_fixed[9]) {
			i3_flt r = (max_p->disc_index - min_p->disc_index);
			P->disc_index += scale * r * i3_random_normal();
			if (P->disc_index<min_p->disc_index) P->disc_index = min_p->disc_index+r*epsilon;
			if (P->disc_index>max_p->disc_index) P->disc_index = max_p->disc_index-r*epsilon;
		}
	}
	if (!strcmp(model->name,"multixray")){
		i3_multixray_parameter_set * P = (i3_multixray_parameter_set*) p;
		i3_multixray_parameter_set * min_p = (i3_multixray_parameter_set*)model->min;
		i3_multixray_parameter_set * max_p = (i3_multixray_parameter_set*)model->max;
		i3_flt epsilon = 1e-6;
		 if (!model->param_fixed[0]) {
			i3_flt r = (max_p->alpha_arcsec - min_p->alpha_arcsec);
			P->alpha_arcsec += scale * r * i3_random_normal();
			if (P->alpha_arcsec<min_p->alpha_arcsec) P->alpha_arcsec = min_p->alpha_arcsec+r*epsilon;
			if (P->alpha_arcsec>max_p->alpha_arcsec) P->alpha_arcsec = max_p->alpha_arcsec-r*epsilon;
		}
		 if (!model->param_fixed[1]) {
			i3_flt r = (max_p->delta_arcsec - min_p->delta_arcsec);
			P->delta_arcsec += scale * r * i3_random_normal();
			if (P->delta_arcsec<min_p->delta_arcsec) P->delta_arcsec = min_p->delta_arcsec+r*epsilon;
			if (P->delta_arcsec>max_p->delta_arcsec) P->delta_arcsec = max_p->delta_arcsec-r*epsilon;
		}
		 if (!model->param_fixed[2]) {
			i3_flt r = (max_p->e1 - min_p->e1);
			P->e1 += scale * r * i3_random_normal();
			if (P->e1<min_p->e1) P->e1 = min_p->e1+r*epsilon;
			if (P->e1>max_p->e1) P->e1 = max_p->e1-r*epsilon;
		}
		 if (!model->param_fixed[3]) {
			i3_flt r = (max_p->e2 - min_p->e2);
			P->e2 += scale * r * i3_random_normal();
			if (P->e2<min_p->e2) P->e2 = min_p->e2+r*epsilon;
			if (P->e2>max_p->e2) P->e2 = max_p->e2-r*epsilon;
		}
		 if (!model->param_fixed[4]) {
			i3_flt r = (max_p->rc - min_p->rc);
			P->rc += scale * r * i3_random_normal();
			if (P->rc<min_p->rc) P->rc = min_p->rc+r*epsilon;
			if (P->rc>max_p->rc) P->rc = max_p->rc-r*epsilon;
		}
		 if (!model->param_fixed[5]) {
			i3_flt r = (max_p->A - min_p->A);
			P->A += scale * r * i3_random_normal();
			if (P->A<min_p->A) P->A = min_p->A+r*epsilon;
			if (P->A>max_p->A) P->A = max_p->A-r*epsilon;
		}
		 if (!model->param_fixed[6]) {
			i3_flt r = (max_p->beta - min_p->beta);
			P->beta += scale * r * i3_random_normal();
			if (P->beta<min_p->beta) P->beta = min_p->beta+r*epsilon;
			if (P->beta>max_p->beta) P->beta = max_p->beta-r*epsilon;
		}
		 if (!model->param_fixed[7]) {
			 for (int ip=0; ip<30; ip++){
				i3_flt r = (max_p->background[ip] - min_p->background[ip]);
				P->background[ip] += scale * r * i3_random_normal();
				if (P->background[ip]<min_p->background[ip]) P->background[ip] = min_p->background[ip]+r*epsilon;
				if (P->background[ip]>max_p->background[ip]) P->background[ip] = max_p->background[ip]-r*epsilon;
			}
		}
	}
}

const char * i3_model_header_line(i3_model * model)
{
		if (!strcmp(model->name,"sersics")){
		return "x0              y0              e1              e2              radius              radius_ratio              bulge_A              disc_A              bulge_index              disc_index              delta_e_bulge              delta_theta_bulge";
	}
	if (!strcmp(model->name,"logsersics")){
		return "x0              y0              e1              e2              radius              radius_ratio              bulge_A              disc_A              bulge_index              disc_index              delta_e_bulge              delta_theta_bulge";
	}
	if (!strcmp(model->name,"quadratic_test")){
		return "x1              x2              x3              x4              x5";
	}
	if (!strcmp(model->name,"multisersics")){
		return "ra_as              dec_as              e1              e2              radius              radius_ratio              bulge_A              disc_A              bulge_index              disc_index";
	}
	if (!strcmp(model->name,"multixray")){
		return "alpha_arcsec              delta_arcsec              e1              e2              rc              A              beta              background";
	}
	return "";
}


void i3_model_pretty_fprintf_parameters(i3_model * model, FILE * f, i3_parameter_set * p){
	char s[i3_max_string_length];
	i3_model_parameter_pretty_string(model,p,s,i3_max_string_length);
	fprintf(f,"%s",s);
}

void i3_model_fprintf_parameters(i3_model * model, FILE * f, i3_parameter_set * p){
	char s[i3_max_string_length];
	i3_model_parameter_string(model,p,s,i3_max_string_length);
	fprintf(f,"%s",s);
}

void i3_model_pretty_printf_parameters(i3_model * model, i3_parameter_set * p){
	i3_model_pretty_fprintf_parameters(model,stdout,p);
}

void i3_model_printf_parameters(i3_model * model, i3_parameter_set * p){
	i3_model_fprintf_parameters(model,stdout,p);
}


int i3_model_parameter_string(i3_model * model, i3_parameter_set * p, char * parameter_string, int string_length){
	if (!strcmp(model->name,"sersics")){
		i3_sersics_parameter_set * mp = (i3_sersics_parameter_set*) p;
		return snprintf(parameter_string,string_length,"% e\t% e\t% e\t% e\t% e\t% e\t% e\t% e\t% e\t% e\t% e\t% e",mp->x0,mp->y0,mp->e1,mp->e2,mp->radius,mp->radius_ratio,mp->bulge_A,mp->disc_A,mp->bulge_index,mp->disc_index,mp->delta_e_bulge,mp->delta_theta_bulge);
	}
	if (!strcmp(model->name,"logsersics")){
		i3_logsersics_parameter_set * mp = (i3_logsersics_parameter_set*) p;
		return snprintf(parameter_string,string_length,"% e\t% e\t% e\t% e\t% e\t% e\t% e\t% e\t% e\t% e\t% e\t% e",mp->x0,mp->y0,mp->e1,mp->e2,mp->radius,mp->radius_ratio,mp->bulge_A,mp->disc_A,mp->bulge_index,mp->disc_index,mp->delta_e_bulge,mp->delta_theta_bulge);
	}
	if (!strcmp(model->name,"quadratic_test")){
		i3_quadratic_test_parameter_set * mp = (i3_quadratic_test_parameter_set*) p;
		return snprintf(parameter_string,string_length,"% e\t% e\t% e\t% e\t%e %e %e ",mp->x1,mp->x2,mp->x3,mp->x4,mp->x5[0], mp->x5[1], mp->x5[2]);
	}
	if (!strcmp(model->name,"multisersics")){
		i3_multisersics_parameter_set * mp = (i3_multisersics_parameter_set*) p;
		return snprintf(parameter_string,string_length,"% e\t% e\t% e\t% e\t% e\t% e\t% e\t% e\t% e\t% e",mp->ra_as,mp->dec_as,mp->e1,mp->e2,mp->radius,mp->radius_ratio,mp->bulge_A,mp->disc_A,mp->bulge_index,mp->disc_index);
	}
	if (!strcmp(model->name,"multixray")){
		i3_multixray_parameter_set * mp = (i3_multixray_parameter_set*) p;
		return snprintf(parameter_string,string_length,"% e\t% e\t% e\t% e\t% e\t% e\t% e\t%e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e ",mp->alpha_arcsec,mp->delta_arcsec,mp->e1,mp->e2,mp->rc,mp->A,mp->beta,mp->background[0], mp->background[1], mp->background[2], mp->background[3], mp->background[4], mp->background[5], mp->background[6], mp->background[7], mp->background[8], mp->background[9], mp->background[10], mp->background[11], mp->background[12], mp->background[13], mp->background[14], mp->background[15], mp->background[16], mp->background[17], mp->background[18], mp->background[19], mp->background[20], mp->background[21], mp->background[22], mp->background[23], mp->background[24], mp->background[25], mp->background[26], mp->background[27], mp->background[28], mp->background[29]);
	}
	return 0;
}
int i3_model_parameter_pretty_string(i3_model * model, i3_parameter_set * p, char * parameter_string, int string_length){
	if (!strcmp(model->name,"sersics")){
		i3_sersics_parameter_set * mp = (i3_sersics_parameter_set*) p;
		return snprintf(parameter_string,string_length,"x0=% e\ny0=% e\ne1=% e\ne2=% e\nradius=% e\nradius_ratio=% e\nbulge_A=% e\ndisc_A=% e\nbulge_index=% e\ndisc_index=% e\ndelta_e_bulge=% e\ndelta_theta_bulge=% e",mp->x0,mp->y0,mp->e1,mp->e2,mp->radius,mp->radius_ratio,mp->bulge_A,mp->disc_A,mp->bulge_index,mp->disc_index,mp->delta_e_bulge,mp->delta_theta_bulge);
	}
	if (!strcmp(model->name,"logsersics")){
		i3_logsersics_parameter_set * mp = (i3_logsersics_parameter_set*) p;
		return snprintf(parameter_string,string_length,"x0=% e\ny0=% e\ne1=% e\ne2=% e\nradius=% e\nradius_ratio=% e\nbulge_A=% e\ndisc_A=% e\nbulge_index=% e\ndisc_index=% e\ndelta_e_bulge=% e\ndelta_theta_bulge=% e",mp->x0,mp->y0,mp->e1,mp->e2,mp->radius,mp->radius_ratio,mp->bulge_A,mp->disc_A,mp->bulge_index,mp->disc_index,mp->delta_e_bulge,mp->delta_theta_bulge);
	}
	if (!strcmp(model->name,"quadratic_test")){
		i3_quadratic_test_parameter_set * mp = (i3_quadratic_test_parameter_set*) p;
		return snprintf(parameter_string,string_length,"x1=% e\nx2=% e\nx3=% e\nx4=% e\nx5=%e, %e, %e",mp->x1,mp->x2,mp->x3,mp->x4,mp->x5[0], mp->x5[1], mp->x5[2]);
	}
	if (!strcmp(model->name,"multisersics")){
		i3_multisersics_parameter_set * mp = (i3_multisersics_parameter_set*) p;
		return snprintf(parameter_string,string_length,"ra_as=% e\ndec_as=% e\ne1=% e\ne2=% e\nradius=% e\nradius_ratio=% e\nbulge_A=% e\ndisc_A=% e\nbulge_index=% e\ndisc_index=% e",mp->ra_as,mp->dec_as,mp->e1,mp->e2,mp->radius,mp->radius_ratio,mp->bulge_A,mp->disc_A,mp->bulge_index,mp->disc_index);
	}
	if (!strcmp(model->name,"multixray")){
		i3_multixray_parameter_set * mp = (i3_multixray_parameter_set*) p;
		return snprintf(parameter_string,string_length,"alpha_arcsec=% e\ndelta_arcsec=% e\ne1=% e\ne2=% e\nrc=% e\nA=% e\nbeta=% e\nbackground=%e, %e, %e, %e, %e, %e, %e, %e, %e, %e, %e, %e, %e, %e, %e, %e, %e, %e, %e, %e, %e, %e, %e, %e, %e, %e, %e, %e, %e, %e",mp->alpha_arcsec,mp->delta_arcsec,mp->e1,mp->e2,mp->rc,mp->A,mp->beta,mp->background[0], mp->background[1], mp->background[2], mp->background[3], mp->background[4], mp->background[5], mp->background[6], mp->background[7], mp->background[8], mp->background[9], mp->background[10], mp->background[11], mp->background[12], mp->background[13], mp->background[14], mp->background[15], mp->background[16], mp->background[17], mp->background[18], mp->background[19], mp->background[20], mp->background[21], mp->background[22], mp->background[23], mp->background[24], mp->background[25], mp->background[26], mp->background[27], mp->background[28], mp->background[29]);
	}
	return 0;
}

void i3_model_map_physical(i3_model * model, i3_parameter_set * input, i3_parameter_set * output, i3_options * options)
{
	i3_model_copy_parameters(model, output, input);
	if (model->map_physical!=NULL) model->map_physical(input, output, options);
}

void i3_model_starting_position(i3_model * model, i3_parameter_set * p, i3_data_set * data_set, i3_options * options){
	if (model->start==NULL) return;
	model->start(p,data_set,options);
}

i3_flt i3_model_posterior(i3_model * model,i3_image * image, i3_parameter_set * p, i3_data_set * data){
	i3_flt prior = i3_model_prior(model,p);
	if (prior==BAD_LIKELIHOOD) return prior;
	i3_flt likelihood=i3_model_likelihood(model, image, p, data);
	if (likelihood==BAD_LIKELIHOOD) return likelihood;
	return prior+likelihood;
}

i3_flt i3_model_likelihood(i3_model * model,i3_image * image, i3_parameter_set * p, i3_data_set * data){
	i3_flt likelihood=model->likelihood(image,p,data);
	return likelihood;
}


i3_flt i3_model_prior(i3_model * model, i3_parameter_set * p){
	i3_flt prior = i3_model_prior_limits(model,p);
	if (prior==BAD_LIKELIHOOD) return prior;  /* Parameters have been excluded because one of them goes outside its range.  No point going further*/
	if (model->prior!=NULL) prior += model->prior(p);
	return prior;
}

i3_flt i3_model_prior_limits(i3_model * model, i3_parameter_set * p){
	if (!strcmp(model->name,"sersics")){
		i3_sersics_parameter_set * P = (i3_sersics_parameter_set*) p;
		i3_sersics_parameter_set * min_p = (i3_sersics_parameter_set*) model->min;
		i3_sersics_parameter_set * max_p = (i3_sersics_parameter_set*) model->max;
		if (P->x0<min_p->x0||P->x0>max_p->x0) return BAD_LIKELIHOOD;
		if (P->y0<min_p->y0||P->y0>max_p->y0) return BAD_LIKELIHOOD;
		if (P->e1<min_p->e1||P->e1>max_p->e1) return BAD_LIKELIHOOD;
		if (P->e2<min_p->e2||P->e2>max_p->e2) return BAD_LIKELIHOOD;
		if (P->radius<min_p->radius||P->radius>max_p->radius) return BAD_LIKELIHOOD;
		if (P->radius_ratio<min_p->radius_ratio||P->radius_ratio>max_p->radius_ratio) return BAD_LIKELIHOOD;
		if (P->bulge_A<min_p->bulge_A||P->bulge_A>max_p->bulge_A) return BAD_LIKELIHOOD;
		if (P->disc_A<min_p->disc_A||P->disc_A>max_p->disc_A) return BAD_LIKELIHOOD;
		if (P->bulge_index<min_p->bulge_index||P->bulge_index>max_p->bulge_index) return BAD_LIKELIHOOD;
		if (P->disc_index<min_p->disc_index||P->disc_index>max_p->disc_index) return BAD_LIKELIHOOD;
		if (P->delta_e_bulge<min_p->delta_e_bulge||P->delta_e_bulge>max_p->delta_e_bulge) return BAD_LIKELIHOOD;
		if (P->delta_theta_bulge<min_p->delta_theta_bulge||P->delta_theta_bulge>max_p->delta_theta_bulge) return BAD_LIKELIHOOD;
	}
	if (!strcmp(model->name,"logsersics")){
		i3_logsersics_parameter_set * P = (i3_logsersics_parameter_set*) p;
		i3_logsersics_parameter_set * min_p = (i3_logsersics_parameter_set*) model->min;
		i3_logsersics_parameter_set * max_p = (i3_logsersics_parameter_set*) model->max;
		if (P->x0<min_p->x0||P->x0>max_p->x0) return BAD_LIKELIHOOD;
		if (P->y0<min_p->y0||P->y0>max_p->y0) return BAD_LIKELIHOOD;
		if (P->e1<min_p->e1||P->e1>max_p->e1) return BAD_LIKELIHOOD;
		if (P->e2<min_p->e2||P->e2>max_p->e2) return BAD_LIKELIHOOD;
		if (P->radius<min_p->radius||P->radius>max_p->radius) return BAD_LIKELIHOOD;
		if (P->radius_ratio<min_p->radius_ratio||P->radius_ratio>max_p->radius_ratio) return BAD_LIKELIHOOD;
		if (P->bulge_A<min_p->bulge_A||P->bulge_A>max_p->bulge_A) return BAD_LIKELIHOOD;
		if (P->disc_A<min_p->disc_A||P->disc_A>max_p->disc_A) return BAD_LIKELIHOOD;
		if (P->bulge_index<min_p->bulge_index||P->bulge_index>max_p->bulge_index) return BAD_LIKELIHOOD;
		if (P->disc_index<min_p->disc_index||P->disc_index>max_p->disc_index) return BAD_LIKELIHOOD;
		if (P->delta_e_bulge<min_p->delta_e_bulge||P->delta_e_bulge>max_p->delta_e_bulge) return BAD_LIKELIHOOD;
		if (P->delta_theta_bulge<min_p->delta_theta_bulge||P->delta_theta_bulge>max_p->delta_theta_bulge) return BAD_LIKELIHOOD;
	}
	if (!strcmp(model->name,"quadratic_test")){
		i3_quadratic_test_parameter_set * P = (i3_quadratic_test_parameter_set*) p;
		i3_quadratic_test_parameter_set * min_p = (i3_quadratic_test_parameter_set*) model->min;
		i3_quadratic_test_parameter_set * max_p = (i3_quadratic_test_parameter_set*) model->max;
		if (P->x1<min_p->x1||P->x1>max_p->x1) return BAD_LIKELIHOOD;
		if (P->x2<min_p->x2||P->x2>max_p->x2) return BAD_LIKELIHOOD;
		if (P->x3<min_p->x3||P->x3>max_p->x3) return BAD_LIKELIHOOD;
		if (P->x4<min_p->x4||P->x4>max_p->x4) return BAD_LIKELIHOOD;
		for (int ip=0; ip<3; ip++) if (P->x5[ip]<min_p->x5[ip]||P->x5[ip]>max_p->x5[ip]) return BAD_LIKELIHOOD;
	}
	if (!strcmp(model->name,"multisersics")){
		i3_multisersics_parameter_set * P = (i3_multisersics_parameter_set*) p;
		i3_multisersics_parameter_set * min_p = (i3_multisersics_parameter_set*) model->min;
		i3_multisersics_parameter_set * max_p = (i3_multisersics_parameter_set*) model->max;
		if (P->ra_as<min_p->ra_as||P->ra_as>max_p->ra_as) return BAD_LIKELIHOOD;
		if (P->dec_as<min_p->dec_as||P->dec_as>max_p->dec_as) return BAD_LIKELIHOOD;
		if (P->e1<min_p->e1||P->e1>max_p->e1) return BAD_LIKELIHOOD;
		if (P->e2<min_p->e2||P->e2>max_p->e2) return BAD_LIKELIHOOD;
		if (P->radius<min_p->radius||P->radius>max_p->radius) return BAD_LIKELIHOOD;
		if (P->radius_ratio<min_p->radius_ratio||P->radius_ratio>max_p->radius_ratio) return BAD_LIKELIHOOD;
		if (P->bulge_A<min_p->bulge_A||P->bulge_A>max_p->bulge_A) return BAD_LIKELIHOOD;
		if (P->disc_A<min_p->disc_A||P->disc_A>max_p->disc_A) return BAD_LIKELIHOOD;
		if (P->bulge_index<min_p->bulge_index||P->bulge_index>max_p->bulge_index) return BAD_LIKELIHOOD;
		if (P->disc_index<min_p->disc_index||P->disc_index>max_p->disc_index) return BAD_LIKELIHOOD;
	}
	if (!strcmp(model->name,"multixray")){
		i3_multixray_parameter_set * P = (i3_multixray_parameter_set*) p;
		i3_multixray_parameter_set * min_p = (i3_multixray_parameter_set*) model->min;
		i3_multixray_parameter_set * max_p = (i3_multixray_parameter_set*) model->max;
		if (P->alpha_arcsec<min_p->alpha_arcsec||P->alpha_arcsec>max_p->alpha_arcsec) return BAD_LIKELIHOOD;
		if (P->delta_arcsec<min_p->delta_arcsec||P->delta_arcsec>max_p->delta_arcsec) return BAD_LIKELIHOOD;
		if (P->e1<min_p->e1||P->e1>max_p->e1) return BAD_LIKELIHOOD;
		if (P->e2<min_p->e2||P->e2>max_p->e2) return BAD_LIKELIHOOD;
		if (P->rc<min_p->rc||P->rc>max_p->rc) return BAD_LIKELIHOOD;
		if (P->A<min_p->A||P->A>max_p->A) return BAD_LIKELIHOOD;
		if (P->beta<min_p->beta||P->beta>max_p->beta) return BAD_LIKELIHOOD;
		for (int ip=0; ip<30; ip++) if (P->background[ip]<min_p->background[ip]||P->background[ip]>max_p->background[ip]) return BAD_LIKELIHOOD;
	}	
	return 0.0;
}

bool i3_model_prior_violations(i3_model * model, i3_parameter_set * p, FILE * output){
	bool violated=false;
	
		if (!strcmp(model->name,"sersics")){
		i3_sersics_parameter_set * P = (i3_sersics_parameter_set*) p;
		i3_sersics_parameter_set * min_p = (i3_sersics_parameter_set*) model->min;
		i3_sersics_parameter_set * max_p = (i3_sersics_parameter_set*) model->max;
		if (P->x0<min_p->x0) {fprintf(output,"Violated minimum on sersics x0: %f < %f\n", P->x0,min_p->x0); violated=true;}
		if (P->x0>max_p->x0) {fprintf(output,"Violated maximum on sersics x0: %f > %f\n", P->x0,max_p->x0); violated=true;}
		if (P->y0<min_p->y0) {fprintf(output,"Violated minimum on sersics y0: %f < %f\n", P->y0,min_p->y0); violated=true;}
		if (P->y0>max_p->y0) {fprintf(output,"Violated maximum on sersics y0: %f > %f\n", P->y0,max_p->y0); violated=true;}
		if (P->e1<min_p->e1) {fprintf(output,"Violated minimum on sersics e1: %f < %f\n", P->e1,min_p->e1); violated=true;}
		if (P->e1>max_p->e1) {fprintf(output,"Violated maximum on sersics e1: %f > %f\n", P->e1,max_p->e1); violated=true;}
		if (P->e2<min_p->e2) {fprintf(output,"Violated minimum on sersics e2: %f < %f\n", P->e2,min_p->e2); violated=true;}
		if (P->e2>max_p->e2) {fprintf(output,"Violated maximum on sersics e2: %f > %f\n", P->e2,max_p->e2); violated=true;}
		if (P->radius<min_p->radius) {fprintf(output,"Violated minimum on sersics radius: %f < %f\n", P->radius,min_p->radius); violated=true;}
		if (P->radius>max_p->radius) {fprintf(output,"Violated maximum on sersics radius: %f > %f\n", P->radius,max_p->radius); violated=true;}
		if (P->radius_ratio<min_p->radius_ratio) {fprintf(output,"Violated minimum on sersics radius_ratio: %f < %f\n", P->radius_ratio,min_p->radius_ratio); violated=true;}
		if (P->radius_ratio>max_p->radius_ratio) {fprintf(output,"Violated maximum on sersics radius_ratio: %f > %f\n", P->radius_ratio,max_p->radius_ratio); violated=true;}
		if (P->bulge_A<min_p->bulge_A) {fprintf(output,"Violated minimum on sersics bulge_A: %f < %f\n", P->bulge_A,min_p->bulge_A); violated=true;}
		if (P->bulge_A>max_p->bulge_A) {fprintf(output,"Violated maximum on sersics bulge_A: %f > %f\n", P->bulge_A,max_p->bulge_A); violated=true;}
		if (P->disc_A<min_p->disc_A) {fprintf(output,"Violated minimum on sersics disc_A: %f < %f\n", P->disc_A,min_p->disc_A); violated=true;}
		if (P->disc_A>max_p->disc_A) {fprintf(output,"Violated maximum on sersics disc_A: %f > %f\n", P->disc_A,max_p->disc_A); violated=true;}
		if (P->bulge_index<min_p->bulge_index) {fprintf(output,"Violated minimum on sersics bulge_index: %f < %f\n", P->bulge_index,min_p->bulge_index); violated=true;}
		if (P->bulge_index>max_p->bulge_index) {fprintf(output,"Violated maximum on sersics bulge_index: %f > %f\n", P->bulge_index,max_p->bulge_index); violated=true;}
		if (P->disc_index<min_p->disc_index) {fprintf(output,"Violated minimum on sersics disc_index: %f < %f\n", P->disc_index,min_p->disc_index); violated=true;}
		if (P->disc_index>max_p->disc_index) {fprintf(output,"Violated maximum on sersics disc_index: %f > %f\n", P->disc_index,max_p->disc_index); violated=true;}
		if (P->delta_e_bulge<min_p->delta_e_bulge) {fprintf(output,"Violated minimum on sersics delta_e_bulge: %f < %f\n", P->delta_e_bulge,min_p->delta_e_bulge); violated=true;}
		if (P->delta_e_bulge>max_p->delta_e_bulge) {fprintf(output,"Violated maximum on sersics delta_e_bulge: %f > %f\n", P->delta_e_bulge,max_p->delta_e_bulge); violated=true;}
		if (P->delta_theta_bulge<min_p->delta_theta_bulge) {fprintf(output,"Violated minimum on sersics delta_theta_bulge: %f < %f\n", P->delta_theta_bulge,min_p->delta_theta_bulge); violated=true;}
		if (P->delta_theta_bulge>max_p->delta_theta_bulge) {fprintf(output,"Violated maximum on sersics delta_theta_bulge: %f > %f\n", P->delta_theta_bulge,max_p->delta_theta_bulge); violated=true;}
	}
	if (!strcmp(model->name,"logsersics")){
		i3_logsersics_parameter_set * P = (i3_logsersics_parameter_set*) p;
		i3_logsersics_parameter_set * min_p = (i3_logsersics_parameter_set*) model->min;
		i3_logsersics_parameter_set * max_p = (i3_logsersics_parameter_set*) model->max;
		if (P->x0<min_p->x0) {fprintf(output,"Violated minimum on logsersics x0: %f < %f\n", P->x0,min_p->x0); violated=true;}
		if (P->x0>max_p->x0) {fprintf(output,"Violated maximum on logsersics x0: %f > %f\n", P->x0,max_p->x0); violated=true;}
		if (P->y0<min_p->y0) {fprintf(output,"Violated minimum on logsersics y0: %f < %f\n", P->y0,min_p->y0); violated=true;}
		if (P->y0>max_p->y0) {fprintf(output,"Violated maximum on logsersics y0: %f > %f\n", P->y0,max_p->y0); violated=true;}
		if (P->e1<min_p->e1) {fprintf(output,"Violated minimum on logsersics e1: %f < %f\n", P->e1,min_p->e1); violated=true;}
		if (P->e1>max_p->e1) {fprintf(output,"Violated maximum on logsersics e1: %f > %f\n", P->e1,max_p->e1); violated=true;}
		if (P->e2<min_p->e2) {fprintf(output,"Violated minimum on logsersics e2: %f < %f\n", P->e2,min_p->e2); violated=true;}
		if (P->e2>max_p->e2) {fprintf(output,"Violated maximum on logsersics e2: %f > %f\n", P->e2,max_p->e2); violated=true;}
		if (P->radius<min_p->radius) {fprintf(output,"Violated minimum on logsersics radius: %f < %f\n", P->radius,min_p->radius); violated=true;}
		if (P->radius>max_p->radius) {fprintf(output,"Violated maximum on logsersics radius: %f > %f\n", P->radius,max_p->radius); violated=true;}
		if (P->radius_ratio<min_p->radius_ratio) {fprintf(output,"Violated minimum on logsersics radius_ratio: %f < %f\n", P->radius_ratio,min_p->radius_ratio); violated=true;}
		if (P->radius_ratio>max_p->radius_ratio) {fprintf(output,"Violated maximum on logsersics radius_ratio: %f > %f\n", P->radius_ratio,max_p->radius_ratio); violated=true;}
		if (P->bulge_A<min_p->bulge_A) {fprintf(output,"Violated minimum on logsersics bulge_A: %f < %f\n", P->bulge_A,min_p->bulge_A); violated=true;}
		if (P->bulge_A>max_p->bulge_A) {fprintf(output,"Violated maximum on logsersics bulge_A: %f > %f\n", P->bulge_A,max_p->bulge_A); violated=true;}
		if (P->disc_A<min_p->disc_A) {fprintf(output,"Violated minimum on logsersics disc_A: %f < %f\n", P->disc_A,min_p->disc_A); violated=true;}
		if (P->disc_A>max_p->disc_A) {fprintf(output,"Violated maximum on logsersics disc_A: %f > %f\n", P->disc_A,max_p->disc_A); violated=true;}
		if (P->bulge_index<min_p->bulge_index) {fprintf(output,"Violated minimum on logsersics bulge_index: %f < %f\n", P->bulge_index,min_p->bulge_index); violated=true;}
		if (P->bulge_index>max_p->bulge_index) {fprintf(output,"Violated maximum on logsersics bulge_index: %f > %f\n", P->bulge_index,max_p->bulge_index); violated=true;}
		if (P->disc_index<min_p->disc_index) {fprintf(output,"Violated minimum on logsersics disc_index: %f < %f\n", P->disc_index,min_p->disc_index); violated=true;}
		if (P->disc_index>max_p->disc_index) {fprintf(output,"Violated maximum on logsersics disc_index: %f > %f\n", P->disc_index,max_p->disc_index); violated=true;}
		if (P->delta_e_bulge<min_p->delta_e_bulge) {fprintf(output,"Violated minimum on logsersics delta_e_bulge: %f < %f\n", P->delta_e_bulge,min_p->delta_e_bulge); violated=true;}
		if (P->delta_e_bulge>max_p->delta_e_bulge) {fprintf(output,"Violated maximum on logsersics delta_e_bulge: %f > %f\n", P->delta_e_bulge,max_p->delta_e_bulge); violated=true;}
		if (P->delta_theta_bulge<min_p->delta_theta_bulge) {fprintf(output,"Violated minimum on logsersics delta_theta_bulge: %f < %f\n", P->delta_theta_bulge,min_p->delta_theta_bulge); violated=true;}
		if (P->delta_theta_bulge>max_p->delta_theta_bulge) {fprintf(output,"Violated maximum on logsersics delta_theta_bulge: %f > %f\n", P->delta_theta_bulge,max_p->delta_theta_bulge); violated=true;}
	}
	if (!strcmp(model->name,"quadratic_test")){
		i3_quadratic_test_parameter_set * P = (i3_quadratic_test_parameter_set*) p;
		i3_quadratic_test_parameter_set * min_p = (i3_quadratic_test_parameter_set*) model->min;
		i3_quadratic_test_parameter_set * max_p = (i3_quadratic_test_parameter_set*) model->max;
		if (P->x1<min_p->x1) {fprintf(output,"Violated minimum on quadratic_test x1: %f < %f\n", P->x1,min_p->x1); violated=true;}
		if (P->x1>max_p->x1) {fprintf(output,"Violated maximum on quadratic_test x1: %f > %f\n", P->x1,max_p->x1); violated=true;}
		if (P->x2<min_p->x2) {fprintf(output,"Violated minimum on quadratic_test x2: %f < %f\n", P->x2,min_p->x2); violated=true;}
		if (P->x2>max_p->x2) {fprintf(output,"Violated maximum on quadratic_test x2: %f > %f\n", P->x2,max_p->x2); violated=true;}
		if (P->x3<min_p->x3) {fprintf(output,"Violated minimum on quadratic_test x3: %f < %f\n", P->x3,min_p->x3); violated=true;}
		if (P->x3>max_p->x3) {fprintf(output,"Violated maximum on quadratic_test x3: %f > %f\n", P->x3,max_p->x3); violated=true;}
		if (P->x4<min_p->x4) {fprintf(output,"Violated minimum on quadratic_test x4: %f < %f\n", P->x4,min_p->x4); violated=true;}
		if (P->x4>max_p->x4) {fprintf(output,"Violated maximum on quadratic_test x4: %f > %f\n", P->x4,max_p->x4); violated=true;}
		for (int ip=0; ip<3; ip++) if (P->x5[ip]<min_p->x5[ip]) {fprintf(output,"Violated minimum on quadratic_test x5[%d]: %f < %f\n", ip, P->x5[ip],min_p->x5[ip]); violated=true;}
		for (int ip=0; ip<3; ip++) if (P->x5[ip]>max_p->x5[ip]) {fprintf(output,"Violated maximum on quadratic_test x5[%d]: %f > %f\n", ip, P->x5[ip],max_p->x5[ip]); violated=true;}
	}
	if (!strcmp(model->name,"multisersics")){
		i3_multisersics_parameter_set * P = (i3_multisersics_parameter_set*) p;
		i3_multisersics_parameter_set * min_p = (i3_multisersics_parameter_set*) model->min;
		i3_multisersics_parameter_set * max_p = (i3_multisersics_parameter_set*) model->max;
		if (P->ra_as<min_p->ra_as) {fprintf(output,"Violated minimum on multisersics ra_as: %f < %f\n", P->ra_as,min_p->ra_as); violated=true;}
		if (P->ra_as>max_p->ra_as) {fprintf(output,"Violated maximum on multisersics ra_as: %f > %f\n", P->ra_as,max_p->ra_as); violated=true;}
		if (P->dec_as<min_p->dec_as) {fprintf(output,"Violated minimum on multisersics dec_as: %f < %f\n", P->dec_as,min_p->dec_as); violated=true;}
		if (P->dec_as>max_p->dec_as) {fprintf(output,"Violated maximum on multisersics dec_as: %f > %f\n", P->dec_as,max_p->dec_as); violated=true;}
		if (P->e1<min_p->e1) {fprintf(output,"Violated minimum on multisersics e1: %f < %f\n", P->e1,min_p->e1); violated=true;}
		if (P->e1>max_p->e1) {fprintf(output,"Violated maximum on multisersics e1: %f > %f\n", P->e1,max_p->e1); violated=true;}
		if (P->e2<min_p->e2) {fprintf(output,"Violated minimum on multisersics e2: %f < %f\n", P->e2,min_p->e2); violated=true;}
		if (P->e2>max_p->e2) {fprintf(output,"Violated maximum on multisersics e2: %f > %f\n", P->e2,max_p->e2); violated=true;}
		if (P->radius<min_p->radius) {fprintf(output,"Violated minimum on multisersics radius: %f < %f\n", P->radius,min_p->radius); violated=true;}
		if (P->radius>max_p->radius) {fprintf(output,"Violated maximum on multisersics radius: %f > %f\n", P->radius,max_p->radius); violated=true;}
		if (P->radius_ratio<min_p->radius_ratio) {fprintf(output,"Violated minimum on multisersics radius_ratio: %f < %f\n", P->radius_ratio,min_p->radius_ratio); violated=true;}
		if (P->radius_ratio>max_p->radius_ratio) {fprintf(output,"Violated maximum on multisersics radius_ratio: %f > %f\n", P->radius_ratio,max_p->radius_ratio); violated=true;}
		if (P->bulge_A<min_p->bulge_A) {fprintf(output,"Violated minimum on multisersics bulge_A: %f < %f\n", P->bulge_A,min_p->bulge_A); violated=true;}
		if (P->bulge_A>max_p->bulge_A) {fprintf(output,"Violated maximum on multisersics bulge_A: %f > %f\n", P->bulge_A,max_p->bulge_A); violated=true;}
		if (P->disc_A<min_p->disc_A) {fprintf(output,"Violated minimum on multisersics disc_A: %f < %f\n", P->disc_A,min_p->disc_A); violated=true;}
		if (P->disc_A>max_p->disc_A) {fprintf(output,"Violated maximum on multisersics disc_A: %f > %f\n", P->disc_A,max_p->disc_A); violated=true;}
		if (P->bulge_index<min_p->bulge_index) {fprintf(output,"Violated minimum on multisersics bulge_index: %f < %f\n", P->bulge_index,min_p->bulge_index); violated=true;}
		if (P->bulge_index>max_p->bulge_index) {fprintf(output,"Violated maximum on multisersics bulge_index: %f > %f\n", P->bulge_index,max_p->bulge_index); violated=true;}
		if (P->disc_index<min_p->disc_index) {fprintf(output,"Violated minimum on multisersics disc_index: %f < %f\n", P->disc_index,min_p->disc_index); violated=true;}
		if (P->disc_index>max_p->disc_index) {fprintf(output,"Violated maximum on multisersics disc_index: %f > %f\n", P->disc_index,max_p->disc_index); violated=true;}
	}
	if (!strcmp(model->name,"multixray")){
		i3_multixray_parameter_set * P = (i3_multixray_parameter_set*) p;
		i3_multixray_parameter_set * min_p = (i3_multixray_parameter_set*) model->min;
		i3_multixray_parameter_set * max_p = (i3_multixray_parameter_set*) model->max;
		if (P->alpha_arcsec<min_p->alpha_arcsec) {fprintf(output,"Violated minimum on multixray alpha_arcsec: %f < %f\n", P->alpha_arcsec,min_p->alpha_arcsec); violated=true;}
		if (P->alpha_arcsec>max_p->alpha_arcsec) {fprintf(output,"Violated maximum on multixray alpha_arcsec: %f > %f\n", P->alpha_arcsec,max_p->alpha_arcsec); violated=true;}
		if (P->delta_arcsec<min_p->delta_arcsec) {fprintf(output,"Violated minimum on multixray delta_arcsec: %f < %f\n", P->delta_arcsec,min_p->delta_arcsec); violated=true;}
		if (P->delta_arcsec>max_p->delta_arcsec) {fprintf(output,"Violated maximum on multixray delta_arcsec: %f > %f\n", P->delta_arcsec,max_p->delta_arcsec); violated=true;}
		if (P->e1<min_p->e1) {fprintf(output,"Violated minimum on multixray e1: %f < %f\n", P->e1,min_p->e1); violated=true;}
		if (P->e1>max_p->e1) {fprintf(output,"Violated maximum on multixray e1: %f > %f\n", P->e1,max_p->e1); violated=true;}
		if (P->e2<min_p->e2) {fprintf(output,"Violated minimum on multixray e2: %f < %f\n", P->e2,min_p->e2); violated=true;}
		if (P->e2>max_p->e2) {fprintf(output,"Violated maximum on multixray e2: %f > %f\n", P->e2,max_p->e2); violated=true;}
		if (P->rc<min_p->rc) {fprintf(output,"Violated minimum on multixray rc: %f < %f\n", P->rc,min_p->rc); violated=true;}
		if (P->rc>max_p->rc) {fprintf(output,"Violated maximum on multixray rc: %f > %f\n", P->rc,max_p->rc); violated=true;}
		if (P->A<min_p->A) {fprintf(output,"Violated minimum on multixray A: %f < %f\n", P->A,min_p->A); violated=true;}
		if (P->A>max_p->A) {fprintf(output,"Violated maximum on multixray A: %f > %f\n", P->A,max_p->A); violated=true;}
		if (P->beta<min_p->beta) {fprintf(output,"Violated minimum on multixray beta: %f < %f\n", P->beta,min_p->beta); violated=true;}
		if (P->beta>max_p->beta) {fprintf(output,"Violated maximum on multixray beta: %f > %f\n", P->beta,max_p->beta); violated=true;}
		for (int ip=0; ip<30; ip++) if (P->background[ip]<min_p->background[ip]) {fprintf(output,"Violated minimum on multixray background[%d]: %f < %f\n", ip, P->background[ip],min_p->background[ip]); violated=true;}
		for (int ip=0; ip<30; ip++) if (P->background[ip]>max_p->background[ip]) {fprintf(output,"Violated maximum on multixray background[%d]: %f > %f\n", ip, P->background[ip],max_p->background[ip]); violated=true;}
	}
	
	return violated;
}

void i3_model_get_image_into(i3_model * model, i3_parameter_set * params, i3_data_set * dataset, i3_image * image_out){
	if (!strcmp(model->name,"sersics")){
		i3_sersics_model_image(params,dataset,image_out);
	}
	if (!strcmp(model->name,"logsersics")){
		i3_logsersics_model_image(params,dataset,image_out);
	}
	if (!strcmp(model->name,"quadratic_test")){
		 printf("image generation function for this model quadratic_test is not implenented.\n");
	}
	if (!strcmp(model->name,"multisersics")){
		 printf("image generation function for this model multisersics is not implenented.\n");
	}
	if (!strcmp(model->name,"multixray")){
		 printf("image generation function for this model multixray is not implenented.\n");
	}
}

void i3_model_read_params_from_str_into(i3_model * model, char *parameter_strings[], i3_parameter_set * params){
	if (!strcmp(model->name,"sersics")){
		i3_sersics_parameter_set * m_params = (i3_sersics_parameter_set*) params;
		m_params->x0 = atof(parameter_strings[0]);
		m_params->y0 = atof(parameter_strings[1]);
		m_params->e1 = atof(parameter_strings[2]);
		m_params->e2 = atof(parameter_strings[3]);
		m_params->radius = atof(parameter_strings[4]);
		m_params->radius_ratio = atof(parameter_strings[5]);
		m_params->bulge_A = atof(parameter_strings[6]);
		m_params->disc_A = atof(parameter_strings[7]);
		m_params->bulge_index = atof(parameter_strings[8]);
		m_params->disc_index = atof(parameter_strings[9]);
		m_params->delta_e_bulge = atof(parameter_strings[10]);
		m_params->delta_theta_bulge = atof(parameter_strings[11]);
	}
	if (!strcmp(model->name,"logsersics")){
		i3_logsersics_parameter_set * m_params = (i3_logsersics_parameter_set*) params;
		m_params->x0 = atof(parameter_strings[0]);
		m_params->y0 = atof(parameter_strings[1]);
		m_params->e1 = atof(parameter_strings[2]);
		m_params->e2 = atof(parameter_strings[3]);
		m_params->radius = atof(parameter_strings[4]);
		m_params->radius_ratio = atof(parameter_strings[5]);
		m_params->bulge_A = atof(parameter_strings[6]);
		m_params->disc_A = atof(parameter_strings[7]);
		m_params->bulge_index = atof(parameter_strings[8]);
		m_params->disc_index = atof(parameter_strings[9]);
		m_params->delta_e_bulge = atof(parameter_strings[10]);
		m_params->delta_theta_bulge = atof(parameter_strings[11]);
	}
	if (!strcmp(model->name,"quadratic_test")){
		i3_quadratic_test_parameter_set * m_params = (i3_quadratic_test_parameter_set*) params;
		m_params->x1 = atof(parameter_strings[0]);
		m_params->x2 = atof(parameter_strings[1]);
		m_params->x3 = atof(parameter_strings[2]);
		m_params->x4 = atof(parameter_strings[3]);
		int status; i3_options_parse_flt_array(m_params->x5,parameter_strings[4],3, &status);
	}
	if (!strcmp(model->name,"multisersics")){
		i3_multisersics_parameter_set * m_params = (i3_multisersics_parameter_set*) params;
		m_params->ra_as = atof(parameter_strings[0]);
		m_params->dec_as = atof(parameter_strings[1]);
		m_params->e1 = atof(parameter_strings[2]);
		m_params->e2 = atof(parameter_strings[3]);
		m_params->radius = atof(parameter_strings[4]);
		m_params->radius_ratio = atof(parameter_strings[5]);
		m_params->bulge_A = atof(parameter_strings[6]);
		m_params->disc_A = atof(parameter_strings[7]);
		m_params->bulge_index = atof(parameter_strings[8]);
		m_params->disc_index = atof(parameter_strings[9]);
	}
	if (!strcmp(model->name,"multixray")){
		i3_multixray_parameter_set * m_params = (i3_multixray_parameter_set*) params;
		m_params->alpha_arcsec = atof(parameter_strings[0]);
		m_params->delta_arcsec = atof(parameter_strings[1]);
		m_params->e1 = atof(parameter_strings[2]);
		m_params->e2 = atof(parameter_strings[3]);
		m_params->rc = atof(parameter_strings[4]);
		m_params->A = atof(parameter_strings[5]);
		m_params->beta = atof(parameter_strings[6]);
		int status; i3_options_parse_flt_array(m_params->background,parameter_strings[7],30, &status);
	}
}

void i3_model_posterior_derivative_approximation(i3_model * model, i3_parameter_set * p, i3_image * model_image, i3_data_set * data_set, i3_parameter_set * pprime, i3_flt * like, i3_flt epsilon){
	if (!strcmp(model->name,"sersics")){
		i3_sersics_parameter_set * P0 = (i3_sersics_parameter_set*) p;
		i3_sersics_parameter_set * Pprime = (i3_sersics_parameter_set*) pprime;
		i3_sersics_parameter_set * P = (i3_sersics_parameter_set*)malloc(model->nbytes);
		i3_flt L0 = i3_model_posterior(model,model_image,P0,data_set);
		i3_model_copy_parameters(model,P,P0);
		P->x0 += epsilon;
		Pprime->x0 = (i3_model_posterior(model,model_image,P,data_set) - L0)/epsilon;
		i3_model_copy_parameters(model,P,P0);
		P->y0 += epsilon;
		Pprime->y0 = (i3_model_posterior(model,model_image,P,data_set) - L0)/epsilon;
		i3_model_copy_parameters(model,P,P0);
		P->e1 += epsilon;
		Pprime->e1 = (i3_model_posterior(model,model_image,P,data_set) - L0)/epsilon;
		i3_model_copy_parameters(model,P,P0);
		P->e2 += epsilon;
		Pprime->e2 = (i3_model_posterior(model,model_image,P,data_set) - L0)/epsilon;
		i3_model_copy_parameters(model,P,P0);
		P->radius += epsilon;
		Pprime->radius = (i3_model_posterior(model,model_image,P,data_set) - L0)/epsilon;
		i3_model_copy_parameters(model,P,P0);
		P->radius_ratio += epsilon;
		Pprime->radius_ratio = (i3_model_posterior(model,model_image,P,data_set) - L0)/epsilon;
		i3_model_copy_parameters(model,P,P0);
		P->bulge_A += epsilon;
		Pprime->bulge_A = (i3_model_posterior(model,model_image,P,data_set) - L0)/epsilon;
		i3_model_copy_parameters(model,P,P0);
		P->disc_A += epsilon;
		Pprime->disc_A = (i3_model_posterior(model,model_image,P,data_set) - L0)/epsilon;
		i3_model_copy_parameters(model,P,P0);
		P->bulge_index += epsilon;
		Pprime->bulge_index = (i3_model_posterior(model,model_image,P,data_set) - L0)/epsilon;
		i3_model_copy_parameters(model,P,P0);
		P->disc_index += epsilon;
		Pprime->disc_index = (i3_model_posterior(model,model_image,P,data_set) - L0)/epsilon;
		i3_model_copy_parameters(model,P,P0);
		P->delta_e_bulge += epsilon;
		Pprime->delta_e_bulge = (i3_model_posterior(model,model_image,P,data_set) - L0)/epsilon;
		i3_model_copy_parameters(model,P,P0);
		P->delta_theta_bulge += epsilon;
		Pprime->delta_theta_bulge = (i3_model_posterior(model,model_image,P,data_set) - L0)/epsilon;
		Pprime->delta_theta_bulge = (i3_model_posterior(model,model_image,P,data_set) - L0)/epsilon;
		free(P);
	*like=L0;
	}
	if (!strcmp(model->name,"logsersics")){
		i3_logsersics_parameter_set * P0 = (i3_logsersics_parameter_set*) p;
		i3_logsersics_parameter_set * Pprime = (i3_logsersics_parameter_set*) pprime;
		i3_logsersics_parameter_set * P = (i3_logsersics_parameter_set*)malloc(model->nbytes);
		i3_flt L0 = i3_model_posterior(model,model_image,P0,data_set);
		i3_model_copy_parameters(model,P,P0);
		P->x0 += epsilon;
		Pprime->x0 = (i3_model_posterior(model,model_image,P,data_set) - L0)/epsilon;
		i3_model_copy_parameters(model,P,P0);
		P->y0 += epsilon;
		Pprime->y0 = (i3_model_posterior(model,model_image,P,data_set) - L0)/epsilon;
		i3_model_copy_parameters(model,P,P0);
		P->e1 += epsilon;
		Pprime->e1 = (i3_model_posterior(model,model_image,P,data_set) - L0)/epsilon;
		i3_model_copy_parameters(model,P,P0);
		P->e2 += epsilon;
		Pprime->e2 = (i3_model_posterior(model,model_image,P,data_set) - L0)/epsilon;
		i3_model_copy_parameters(model,P,P0);
		P->radius += epsilon;
		Pprime->radius = (i3_model_posterior(model,model_image,P,data_set) - L0)/epsilon;
		i3_model_copy_parameters(model,P,P0);
		P->radius_ratio += epsilon;
		Pprime->radius_ratio = (i3_model_posterior(model,model_image,P,data_set) - L0)/epsilon;
		i3_model_copy_parameters(model,P,P0);
		P->bulge_A += epsilon;
		Pprime->bulge_A = (i3_model_posterior(model,model_image,P,data_set) - L0)/epsilon;
		i3_model_copy_parameters(model,P,P0);
		P->disc_A += epsilon;
		Pprime->disc_A = (i3_model_posterior(model,model_image,P,data_set) - L0)/epsilon;
		i3_model_copy_parameters(model,P,P0);
		P->bulge_index += epsilon;
		Pprime->bulge_index = (i3_model_posterior(model,model_image,P,data_set) - L0)/epsilon;
		i3_model_copy_parameters(model,P,P0);
		P->disc_index += epsilon;
		Pprime->disc_index = (i3_model_posterior(model,model_image,P,data_set) - L0)/epsilon;
		i3_model_copy_parameters(model,P,P0);
		P->delta_e_bulge += epsilon;
		Pprime->delta_e_bulge = (i3_model_posterior(model,model_image,P,data_set) - L0)/epsilon;
		i3_model_copy_parameters(model,P,P0);
		P->delta_theta_bulge += epsilon;
		Pprime->delta_theta_bulge = (i3_model_posterior(model,model_image,P,data_set) - L0)/epsilon;
		Pprime->delta_theta_bulge = (i3_model_posterior(model,model_image,P,data_set) - L0)/epsilon;
		free(P);
	*like=L0;
	}
	if (!strcmp(model->name,"quadratic_test")){
		i3_quadratic_test_parameter_set * P0 = (i3_quadratic_test_parameter_set*) p;
		i3_quadratic_test_parameter_set * Pprime = (i3_quadratic_test_parameter_set*) pprime;
		i3_quadratic_test_parameter_set * P = (i3_quadratic_test_parameter_set*)malloc(model->nbytes);
		i3_flt L0 = i3_model_posterior(model,model_image,P0,data_set);
		i3_model_copy_parameters(model,P,P0);
		P->x1 += epsilon;
		Pprime->x1 = (i3_model_posterior(model,model_image,P,data_set) - L0)/epsilon;
		i3_model_copy_parameters(model,P,P0);
		P->x2 += epsilon;
		Pprime->x2 = (i3_model_posterior(model,model_image,P,data_set) - L0)/epsilon;
		i3_model_copy_parameters(model,P,P0);
		P->x3 += epsilon;
		Pprime->x3 = (i3_model_posterior(model,model_image,P,data_set) - L0)/epsilon;
		i3_model_copy_parameters(model,P,P0);
		P->x4 += epsilon;
		Pprime->x4 = (i3_model_posterior(model,model_image,P,data_set) - L0)/epsilon;
		Pprime->x4 = (i3_model_posterior(model,model_image,P,data_set) - L0)/epsilon;
		free(P);
	*like=L0;
	}
	if (!strcmp(model->name,"multisersics")){
		i3_multisersics_parameter_set * P0 = (i3_multisersics_parameter_set*) p;
		i3_multisersics_parameter_set * Pprime = (i3_multisersics_parameter_set*) pprime;
		i3_multisersics_parameter_set * P = (i3_multisersics_parameter_set*)malloc(model->nbytes);
		i3_flt L0 = i3_model_posterior(model,model_image,P0,data_set);
		i3_model_copy_parameters(model,P,P0);
		P->ra_as += epsilon;
		Pprime->ra_as = (i3_model_posterior(model,model_image,P,data_set) - L0)/epsilon;
		i3_model_copy_parameters(model,P,P0);
		P->dec_as += epsilon;
		Pprime->dec_as = (i3_model_posterior(model,model_image,P,data_set) - L0)/epsilon;
		i3_model_copy_parameters(model,P,P0);
		P->e1 += epsilon;
		Pprime->e1 = (i3_model_posterior(model,model_image,P,data_set) - L0)/epsilon;
		i3_model_copy_parameters(model,P,P0);
		P->e2 += epsilon;
		Pprime->e2 = (i3_model_posterior(model,model_image,P,data_set) - L0)/epsilon;
		i3_model_copy_parameters(model,P,P0);
		P->radius += epsilon;
		Pprime->radius = (i3_model_posterior(model,model_image,P,data_set) - L0)/epsilon;
		i3_model_copy_parameters(model,P,P0);
		P->radius_ratio += epsilon;
		Pprime->radius_ratio = (i3_model_posterior(model,model_image,P,data_set) - L0)/epsilon;
		i3_model_copy_parameters(model,P,P0);
		P->bulge_A += epsilon;
		Pprime->bulge_A = (i3_model_posterior(model,model_image,P,data_set) - L0)/epsilon;
		i3_model_copy_parameters(model,P,P0);
		P->disc_A += epsilon;
		Pprime->disc_A = (i3_model_posterior(model,model_image,P,data_set) - L0)/epsilon;
		i3_model_copy_parameters(model,P,P0);
		P->bulge_index += epsilon;
		Pprime->bulge_index = (i3_model_posterior(model,model_image,P,data_set) - L0)/epsilon;
		i3_model_copy_parameters(model,P,P0);
		P->disc_index += epsilon;
		Pprime->disc_index = (i3_model_posterior(model,model_image,P,data_set) - L0)/epsilon;
		Pprime->disc_index = (i3_model_posterior(model,model_image,P,data_set) - L0)/epsilon;
		free(P);
	*like=L0;
	}
	if (!strcmp(model->name,"multixray")){
		i3_multixray_parameter_set * P0 = (i3_multixray_parameter_set*) p;
		i3_multixray_parameter_set * Pprime = (i3_multixray_parameter_set*) pprime;
		i3_multixray_parameter_set * P = (i3_multixray_parameter_set*)malloc(model->nbytes);
		i3_flt L0 = i3_model_posterior(model,model_image,P0,data_set);
		i3_model_copy_parameters(model,P,P0);
		P->alpha_arcsec += epsilon;
		Pprime->alpha_arcsec = (i3_model_posterior(model,model_image,P,data_set) - L0)/epsilon;
		i3_model_copy_parameters(model,P,P0);
		P->delta_arcsec += epsilon;
		Pprime->delta_arcsec = (i3_model_posterior(model,model_image,P,data_set) - L0)/epsilon;
		i3_model_copy_parameters(model,P,P0);
		P->e1 += epsilon;
		Pprime->e1 = (i3_model_posterior(model,model_image,P,data_set) - L0)/epsilon;
		i3_model_copy_parameters(model,P,P0);
		P->e2 += epsilon;
		Pprime->e2 = (i3_model_posterior(model,model_image,P,data_set) - L0)/epsilon;
		i3_model_copy_parameters(model,P,P0);
		P->rc += epsilon;
		Pprime->rc = (i3_model_posterior(model,model_image,P,data_set) - L0)/epsilon;
		i3_model_copy_parameters(model,P,P0);
		P->A += epsilon;
		Pprime->A = (i3_model_posterior(model,model_image,P,data_set) - L0)/epsilon;
		i3_model_copy_parameters(model,P,P0);
		P->beta += epsilon;
		Pprime->beta = (i3_model_posterior(model,model_image,P,data_set) - L0)/epsilon;
		Pprime->beta = (i3_model_posterior(model,model_image,P,data_set) - L0)/epsilon;
		free(P);
	*like=L0;
	}
}

void i3_model_extract_varied_nonzero_width_parameters(i3_model * model, i3_parameter_set * x, i3_parameter_set * width, i3_flt * output){
        int j=0;
        int offset; 
        double value;
        for(int i=0;i<model->nparam;i++){
                if (model->param_fixed[i]) continue;
                if (!model->param_type[i]==i3_parameter_type_flt) I3_FATAL("Tried to use i3_model_extract_varied_nonzero_width_parameters on non-float parameter.",I3_MATH_ERROR);
                offset = model->byte_offsets[i];  /* The memory position of the parameter in the parameter object. */
                if (0.0==(i3_flt)(*((i3_flt*)(((char*)width)+offset)))) continue; /* Check if width is zero */
                value = (i3_flt)(*((i3_flt*)(((char*)x)+offset)));  /* The memory position of the parameter in the parameter object. */
                output[j++]=value;
        }
}

void i3_model_extract_varied_parameters(i3_model * model, i3_parameter_set * x, i3_flt * output){
        int j=0;
        int offset; 
        double value;
        for(int i=0;i<model->nparam;i++){
                if (model->param_fixed[i]) continue;
                offset = model->byte_offsets[i];  /* The memory position of the parameter in the parameter object. */
                if (!model->param_type[i]==i3_parameter_type_flt) I3_FATAL("Tried to use extract_varied_parmameters on non-float parameter.",I3_MATH_ERROR);
                value = (i3_flt)(*((i3_flt*)(((char*)x)+offset)));  /* The memory position of the parameter in the parameter object. */
                output[j++]=value;
        }
}



void i3_model_input_varied_nonzero_width_parameters(i3_model * model, i3_flt * input, i3_parameter_set * x, i3_parameter_set * width, i3_parameter_set * base_x){
        i3_model_copy_parameters(model,x,base_x);
        
        int j=0;
        for(int i=0;i<model->nparam;i++){
                if (model->param_fixed[i]) continue;
                if (!model->param_type[i]==i3_parameter_type_flt) I3_FATAL("Tried to use i3_model_input_varied_nonzero_width_parameters on non-float parameter.",I3_MATH_ERROR);
                int offset = model->byte_offsets[i];
                if (0.0==(i3_flt)(*((i3_flt*)(((char*)width)+offset)))) continue; /* Check if width is zero */
                *((i3_flt*)(((char*)x)+offset)) = input[j++];
        }
        
}

void i3_model_input_varied_parameters(i3_model * model, i3_flt * input, i3_parameter_set * x, i3_parameter_set * base_x){
        i3_model_copy_parameters(model,x,base_x);
        
        int j=0;
        for(int i=0;i<model->nparam;i++){
                if (model->param_fixed[i]) continue;
                if (!model->param_type[i]==i3_parameter_type_flt) I3_FATAL("Tried to use i3_model_input_varied_parameters on non-float parameter.",I3_MATH_ERROR);
                int offset = model->byte_offsets[i];
                *((i3_flt*)(((char*)x)+offset)) = input[j++];
        }
        
}