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

split files

Browse Source
This commit is contained in:
Yingjie Wang
2025-11-05 12:59:14 -05:00
parent 0330b6b363
commit 8eeb53b9d7
3 changed files with 96 additions and 93 deletions
Download Patch File Download Diff File Expand all files Collapse all files

94
src/background.c
View File

@@ -1,100 +1,8 @@
#include <gsl/gsl_matrix_double.h>
#include <stdint.h>
#include <gsl/gsl_odeiv2.h>
#include <gsl/gsl_matrix.h>
#include <spng.h>
#include "common.h"
// ODE
// du/dphi = v
// dv/dphi = 3Mu^2 - u
int func(double phi, const double y[], double f[], void *params){
f[0] = y[1];
f[1] = -y[0] + 3.0*M*y[0]*y[0];
return GSL_SUCCESS;
}
int jac(double phi, const double y[], double *dfdy, double dfdphi[], void *params){
gsl_matrix_view dfdy_mat = gsl_matrix_view_array(dfdy, 2, 2);
gsl_matrix *m = &dfdy_mat.matrix;
gsl_matrix_set(m, 0, 0, 0.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, 1, 0.0);
dfdphi[0] = 0.0;
dfdphi[0] = 0.0;
return GSL_SUCCESS;
}
double find_root(double *y, double phi){
//printf("phi: %g, y:{%g, %g}\n", phi, y[0], y[1]);
return phi - y[0]/y[1];
}
double find_root_4ord(double u1, double u2, double v1, double v2, double phi1, double phi2){
double Du = u1 - u2;
double Dphi = phi2 - phi1;
return -(u2*v1*Du*Du*Du - v1*v2*Du*Du*Du*phi2 + u2*u2*Du*(2*v1*Du + v2*Du + 3*v1*v2*Dphi) + u2*u2*u2*(v2*Du + v1*(Du + 2*v2*Dphi)))/(v1*v2*Du*Du*Du);
}
inline double find_root_5ord(double u1, double u2, double v1, double v2, double phi1, double phi2){
return (-(phi1*pow(u2,5)*v1*v2) +
u1*pow(u2,4)*(u2 - phi1*v1)*
v2 + pow(u1,2)*pow(u2,3)*
(u2 + 8*phi1*v1)*v2 +
pow(u1,4)*u2*v1*
(-u2 + phi2*v2) +
2*pow(u1,3)*pow(u2,2)*
(u2*(v1 - v2) - 4*phi2*v1*v2)
+ pow(u1,5)*
(-(u2*v1) + phi2*v1*v2))/
(pow(u1 - u2,3)*
(pow(u1,2) + 4*u1*u2 +
pow(u2,2))*v1*v2);
}
// \chi(b) is the totol deflection angle with parameter b
double chi(double b){
double x;
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);
double u0 = 1.0/R0;
double y[2] = {
u0,
sqrt(1.0/(b*b) - u0*u0 + 2.0*M*u0*u0*u0)
};
double phi = 0;
double phi_max = 8*PI;
double y_previous[2] = {0.0, 0.0};
double phi_previous=phi;
while (y[0] > 0 && phi < phi_max){
double phi_next = phi + 0.1;
phi_previous = phi;
y_previous[0] = y[0];
y_previous[1] = y[1];
int status = gsl_odeiv2_driver_apply(d, &phi, phi_next, y);
if (status != GSL_SUCCESS){
printf("Error: GSL driver returns %d\n", status);
break;
}
//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("%.16g\n",find_root_4ord(y_previous[0], y[0], y_previous[1], y[1], phi_previous, phi));
//printf("%.16g\n",find_root_5ord(y_previous[0], y[0], y_previous[1], y[1], phi_previous, phi));
//return (find_root(y_previous, phi_previous) + find_root(y, phi))/2.0;
gsl_odeiv2_driver_free(d);
if (phi >= phi_max){
return phi_max;
}
return find_root_4ord(y_previous[0], y[0], y_previous[1], y[1], phi_previous, phi);
}
#include "init.h"
// (x,y) on screen to far away angle
double phi(double x, double y){

94
src/init.c Normal file
View File

@@ -0,0 +1,94 @@
#include <gsl/gsl_matrix_double.h>
#include <gsl/gsl_odeiv2.h>
#include <gsl/gsl_matrix.h>
#include "common.h"
// ODE
// du/dphi = v
// dv/dphi = 3Mu^2 - u
static int func(double phi, const double y[], double f[], void *params){
f[0] = y[1];
f[1] = -y[0] + 3.0*M*y[0]*y[0];
return GSL_SUCCESS;
}
static int jac(double phi, const double y[], double *dfdy, double dfdphi[], void *params){
gsl_matrix_view dfdy_mat = gsl_matrix_view_array(dfdy, 2, 2);
gsl_matrix *m = &dfdy_mat.matrix;
gsl_matrix_set(m, 0, 0, 0.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, 1, 0.0);
dfdphi[0] = 0.0;
dfdphi[0] = 0.0;
return GSL_SUCCESS;
}
double find_root(double *y, double phi){
//printf("phi: %g, y:{%g, %g}\n", phi, y[0], y[1]);
return phi - y[0]/y[1];
}
double find_root_4ord(double u1, double u2, double v1, double v2, double phi1, double phi2){
double Du = u1 - u2;
double Dphi = phi2 - phi1;
return -(u2*v1*Du*Du*Du - v1*v2*Du*Du*Du*phi2 + u2*u2*Du*(2*v1*Du + v2*Du + 3*v1*v2*Dphi) + u2*u2*u2*(v2*Du + v1*(Du + 2*v2*Dphi)))/(v1*v2*Du*Du*Du);
}
static inline double find_root_5ord(double u1, double u2, double v1, double v2, double phi1, double phi2){
return (-(phi1*pow(u2,5)*v1*v2) +
u1*pow(u2,4)*(u2 - phi1*v1)*
v2 + pow(u1,2)*pow(u2,3)*
(u2 + 8*phi1*v1)*v2 +
pow(u1,4)*u2*v1*
(-u2 + phi2*v2) +
2*pow(u1,3)*pow(u2,2)*
(u2*(v1 - v2) - 4*phi2*v1*v2)
+ pow(u1,5)*
(-(u2*v1) + phi2*v1*v2))/
(pow(u1 - u2,3)*
(pow(u1,2) + 4*u1*u2 +
pow(u2,2))*v1*v2);
}
// \chi(b) is the totol deflection angle with parameter b
double chi(double b){
double x;
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);
double u0 = 1.0/R0;
double y[2] = {
u0,
sqrt(1.0/(b*b) - u0*u0 + 2.0*M*u0*u0*u0)
};
double phi = 0;
double phi_max = 8*PI;
double y_previous[2] = {0.0, 0.0};
double phi_previous=phi;
while (y[0] > 0 && phi < phi_max){
double phi_next = phi + 0.1;
phi_previous = phi;
y_previous[0] = y[0];
y_previous[1] = y[1];
int status = gsl_odeiv2_driver_apply(d, &phi, phi_next, y);
if (status != GSL_SUCCESS){
printf("Error: GSL driver returns %d\n", status);
break;
}
//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("%.16g\n",find_root_4ord(y_previous[0], y[0], y_previous[1], y[1], phi_previous, phi));
//printf("%.16g\n",find_root_5ord(y_previous[0], y[0], y_previous[1], y[1], phi_previous, phi));
//return (find_root(y_previous, phi_previous) + find_root(y, phi))/2.0;
gsl_odeiv2_driver_free(d);
if (phi >= phi_max){
return phi_max;
}
return find_root_5ord(y_previous[0], y[0], y_previous[1], y[1], phi_previous, phi);
}

1
src/init.h Normal file
View File

@@ -0,0 +1 @@
double chi(double b);