fix: if cotpsi is larger than cotpsi_max, color is black

This commit is contained in:
2025-12-07 05:36:55 -05:00
parent 1e1c38da4e
commit 1819d58778
3 changed files with 4 additions and 1 deletions

View File

@@ -9,6 +9,7 @@
#define R0 (15*Rs)
#define f(r) ((1.0-Rs/((double)r)))
#define bmin (1.5*sqrt(3.0)+5e-10)
#define cotpsi_max (sqrt(R0*R0/(bmin*bmin*f(R0)) - 1))
#define THETAERROR 100000
typedef struct {

View File

@@ -145,7 +145,6 @@ int init(double rmax, double rela_err_limit, int *size, double **x, double **y){
double cotpsi_min = 1/rmax;
sample.x[0] = cotpsi_min;
sample.y[0] = chi(cotpsi_min);
double cotpsi_max = sqrt(R0*R0/(bmin*bmin*f(R0)) - 1);
double chi_cotpsimax = chi(cotpsi_max);
refine_interval(cotpsi_min, cotpsi_max, sample.y[0], chi_cotpsimax, rela_err_limit, &sample, 1);

View File

@@ -19,6 +19,7 @@ inline double xy_to_b(double x, double y) {
double integrand_r(double *xy, size_t dim, void *integrand_params_void) {
double cotpsi = 1.0/hypot(xy[0], xy[1]);
if (cotpsi >= cotpsi_max) return 0;
Integrand_params *integrand_params = (Integrand_params *)integrand_params_void;
double theta = gsl_spline_eval(integrand_params->spline_data->spline, cotpsi, integrand_params->spline_data->acc);
double phi = atan2(xy[1], xy[0]);
@@ -30,6 +31,7 @@ double integrand_r(double *xy, size_t dim, void *integrand_params_void) {
double integrand_g(double *xy, size_t dim, void *integrand_params_void) {
double cotpsi = 1.0/hypot(xy[0], xy[1]);
if (cotpsi >= cotpsi_max) return 0;
Integrand_params *integrand_params = (Integrand_params *)integrand_params_void;
double theta = gsl_spline_eval(integrand_params->spline_data->spline, cotpsi, integrand_params->spline_data->acc);
double phi = atan2(xy[1], xy[0]);
@@ -41,6 +43,7 @@ double integrand_g(double *xy, size_t dim, void *integrand_params_void) {
double integrand_b(double *xy, size_t dim, void *integrand_params_void) {
double cotpsi = 1.0/hypot(xy[0], xy[1]);
if (cotpsi >= cotpsi_max) return 0;
Integrand_params *integrand_params = (Integrand_params *)integrand_params_void;
double theta = gsl_spline_eval(integrand_params->spline_data->spline, cotpsi, integrand_params->spline_data->acc);
double phi = atan2(xy[1], xy[0]);