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

fix bug and cleanup

Browse Source
This commit is contained in:
Yingjie Wang
2025-11-10 15:14:21 -05:00
parent d510f37b2e
commit b9954ff8fa
1 changed files with 18 additions and 14 deletions
Download Patch File Download Diff File Expand all files Collapse all files

32
src/init.c
View File

@@ -1,8 +1,11 @@
#include <gsl/gsl_matrix_double.h> #include <gsl/gsl_matrix_double.h>
#include <gsl/gsl_odeiv2.h> #include <gsl/gsl_odeiv2.h>
#include <gsl/gsl_matrix.h> #include <gsl/gsl_matrix.h>
#include <stdio.h>
#include "common.h" #include "common.h"
#define MAX_DEPTH 15
typedef struct { typedef struct {
int capacity; int capacity;
int size; int size;
@@ -24,11 +27,11 @@ static int jac(double phi, const double y[], double *dfdy, double dfdphi[], void
gsl_matrix *m = &dfdy_mat.matrix; gsl_matrix *m = &dfdy_mat.matrix;
gsl_matrix_set(m, 0, 0, 0.0); gsl_matrix_set(m, 0, 0, 0.0);
gsl_matrix_set(m, 0, 1, 1.0); gsl_matrix_set(m, 0, 1, 1.0);
gsl_matrix_set(m, 1, 0, -3*M*y[0]-1.0); gsl_matrix_set(m, 1, 0, 63*M*y[0]-1.0);
gsl_matrix_set(m, 1, 1, 0.0); gsl_matrix_set(m, 1, 1, 0.0);
dfdphi[0] = 0.0; dfdphi[0] = 0.0;
dfdphi[0] = 0.0; dfdphi[1] = 0.0;
return GSL_SUCCESS; return GSL_SUCCESS;
} }
@@ -61,8 +64,7 @@ static inline double find_root_5ord(double u1, double u2, double v1, double v2,
// \chi(b) is the totol deflection angle with parameter b // \chi(b) is the totol deflection angle with parameter b
double chi(double b){ double chi(double b){
double x; gsl_odeiv2_system sys = {func, jac, 2, NULL};
gsl_odeiv2_system sys = {func, jac, 2, &x};
gsl_odeiv2_driver *d = gsl_odeiv2_driver_alloc_y_new(&sys, gsl_odeiv2_step_rk8pd, 1e-3, 1e-12, 1e-10); 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 u0 = 1.0/R0;
@@ -85,7 +87,7 @@ double chi(double b){
printf("Error: GSL driver returns %d\n", status); printf("Error: GSL driver returns %d\n", status);
break; break;
} }
printf("phi: %g, y[0]: %g, y[1]: %g\n",phi, y[0], y[1]); //printf("phi: %g, y[0]: %g, y[1]: %g\n",phi, y[0], y[1]);
} }
//printf("%g, %g \n",find_root(y_previous, phi_previous),find_root(y, phi)); //printf("%g, %g \n",find_root(y_previous, phi_previous),find_root(y, phi));
@@ -106,25 +108,26 @@ void sample_push(SampleData *sample, double x, double y){
sample->x = realloc(sample->x, sizeof(double) * newcap); sample->x = realloc(sample->x, sizeof(double) * newcap);
sample->y = realloc(sample->y, sizeof(double) * newcap); sample->y = realloc(sample->y, sizeof(double) * newcap);
sample->capacity = newcap; sample->capacity = newcap;
//printf("enlarge sample capacity: %d\n", newcap);
} }
sample->x[sample->size] = x; sample->x[sample->size] = x;
sample->y[sample->size] = y; sample->y[sample->size] = y;
sample->size += 1; sample->size += 1;
} }
static void refine_interval(double xl, double xr, double yl, double yr, double rela_err_limit, SampleData *sample){ static void refine_interval(double xl, double xr, double yl, double yr, double rela_err_limit, SampleData *sample, int depth){
printf("calculating interval (%g, %g)\n", xl, xr); //printf("calculating interval (%.16g, %.16g)\n", xl, xr);
double xm = (xl + xr)/2; double xm = (xl + xr)/2;
double ym = chi(xm); double ym = chi(xm);
int need_refine = (fabs(ym-(yl+yr)/2) >= fabs(rela_err_limit * ym)); int need_refine = (fabs(ym-(yl+yr)/2) >= fabs(rela_err_limit * ym)) && depth <= MAX_DEPTH;
if (need_refine) { if (need_refine) {
printf(" "); //printf(" ");
refine_interval(xl, xm, yl, ym, rela_err_limit, sample); refine_interval(xl, xm, yl, ym, rela_err_limit, sample, depth + 1);
} }
sample_push(sample, xm, ym); sample_push(sample, xm, ym);
if (need_refine) { if (need_refine) {
printf(" "); //printf(" ");
refine_interval(xm, xr, ym, yr, rela_err_limit, sample); refine_interval(xm, xr, ym, yr, rela_err_limit, sample, depth + 1);
} }
return ; return ;
} }
@@ -139,14 +142,15 @@ int init(double bmax, double rela_err_limit){
sample.y[0] = chi(bmin); sample.y[0] = chi(bmin);
double chi_bmax = chi(bmax); double chi_bmax = chi(bmax);
refine_interval(bmin, bmax, sample.y[0], chi_bmax, rela_err_limit, &sample); refine_interval(bmin, bmax, sample.y[0], chi_bmax, rela_err_limit, &sample, 1);
sample_push(&sample, bmax, chi_bmax); sample_push(&sample, bmax, chi_bmax);
//sort //sort
//output //output
printf("Total sample points number: %d\n", sample.size); //printf("Total sample points capacity: %d\n", sample.capacity);
//printf("Total sample points number: %d\n", sample.size);
for (int i = 0; i < sample.size; i++) { for (int i = 0; i < sample.size; i++) {
printf("%.16g %.16g\n", sample.x[i], sample.y[i]); printf("%.16g %.16g\n", sample.x[i], sample.y[i]);
} }