diff --git a/src/init.c b/src/init.c index 5514ce3..83bf260 100644 --- a/src/init.c +++ b/src/init.c @@ -63,15 +63,16 @@ static inline double find_root_5ord(double u1, double u2, double v1, double v2, pow(u2,2))*v1*v2); } -// \chi(b) is the totol deflection angle with parameter b -double chi(double b){ +// \chi(cotpsi) is the totol deflection angle with parameter cotpsi=-dr/(r √f dphi) +double chi(double cotpsi){ gsl_odeiv2_system sys = {func, jac, 2, NULL}; gsl_odeiv2_driver *d = gsl_odeiv2_driver_alloc_y_new(&sys, gsl_odeiv2_step_rk8pd, 1e-3, 1e-12, 1e-10); double u0 = 1.0/R0; double y[2] = { u0, - sqrt(1.0/(b*b) - u0*u0 + 2.0*M*u0*u0*u0) + //sqrt(1.0/(b*b) - u0*u0 + 2.0*M*u0*u0*u0) + u0*sqrt(f(R0))*cotpsi }; double phi = 0; double phi_max = 8*PI; @@ -133,7 +134,7 @@ static void refine_interval(double xl, double xr, double yl, double yr, double r return ; } -int init(double bmax, double rela_err_limit, int *size, double **x, double **y){ +int init(double rmax, double rela_err_limit, int *size, double **x, double **y){ if (*x) free(*x); if (*y) free(*y); SampleData sample; @@ -141,13 +142,15 @@ int init(double bmax, double rela_err_limit, int *size, double **x, double **y){ sample.size=1; sample.x = malloc(sizeof(double) * 1); sample.y = malloc(sizeof(double) * 1); - sample.x[0] = bmin; - sample.y[0] = chi(bmin); - double chi_bmax = chi(bmax); + 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(bmin, bmax, sample.y[0], chi_bmax, rela_err_limit, &sample, 1); + refine_interval(cotpsi_min, cotpsi_max, sample.y[0], chi_cotpsimax, rela_err_limit, &sample, 1); - sample_push(&sample, bmax, chi_bmax); + sample_push(&sample, cotpsi_max, chi_cotpsimax); //sort diff --git a/src/init.h b/src/init.h index e6ee675..ae4a1dd 100644 --- a/src/init.h +++ b/src/init.h @@ -1,2 +1,2 @@ -double chi(double b); -int init(double bmax, double rela_err_limit, int *size, double **x, double **y); +double chi(double cotpsi); +int init(double rmax, double rela_err_limit, int *size, double **x, double **y);