wyj/blackhole-ray-trace
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

update: use chi(cotpsi) instead of chi(b)

Browse Source
This commit is contained in:
Yingjie Wang
2025-12-03 15:55:52 -05:00
parent 2b4b36963d
commit 455f15d309
2 changed files with 14 additions and 11 deletions
Download Patch File Download Diff File Expand all files Collapse all files

21
src/init.c
View File

@@ -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

4
src/init.h
View File

@@ -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);