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

update: multiple improvement

Browse Source
This commit is contained in:
Yingjie Wang
2025-12-12 13:47:23 -05:00
parent de6d3bc04a
commit 4be6add815
4 changed files with 34 additions and 33 deletions
Download Patch File Download Diff File Expand all files Collapse all files

6
src/common.h
View File

@@ -5,20 +5,22 @@
#include <gsl/gsl_errno.h> #include <gsl/gsl_errno.h>
#include <math.h> #include <math.h>
#define POW2(x) (x*x) #define POW2(x) ((x)*(x))
#define FLAT 0 #define FLAT 0
#define PI 3.1415926535897932384626433832795028841971693993751058 #define PI 3.1415926535897932384626433832795028841971693993751058
#define Rs 1 #define Rs 1
#define M (FLAT ? 0 : 0.5*Rs) #define M (FLAT ? 0 : 0.5*Rs)
#define R0 (15*Rs) #define R0 (15*Rs)
#define f(r) ((1.0-(2*M)/((double)r))) #define f(r) ((1.0-(2*M)/((double)(r))))
#define bmin (FLAT ? 1e-10 : 1.5*sqrt(3.0)+5e-10) #define bmin (FLAT ? 1e-10 : 1.5*sqrt(3.0)+5e-10)
#define cotpsi_max (sqrt(R0*R0/(bmin*bmin*f(R0)) - 1)) #define cotpsi_max (sqrt(R0*R0/(bmin*bmin*f(R0)) - 1))
#define tanpsi_min (1.0/cotpsi_max)
#define THETAERROR 100000 #define THETAERROR 100000
#define color_index(i,j,c) (j*W*3+i*3+c) #define color_index(i,j,c) (j*W*3+i*3+c)
#define cutperc (FLAT ? 1.0 : 0.99) #define cutperc (FLAT ? 1.0 : 0.99)
#define SCALE 1
typedef struct { typedef struct {
int W; int W;

17
src/init.c
View File

@@ -63,8 +63,8 @@ static inline double find_root_5ord(double u1, double u2, double v1, double v2,
pow(u2,2))*v1*v2); pow(u2,2))*v1*v2);
} }
// \chi(cotpsi) is the totol deflection angle with parameter cotpsi=-dr/(r √f dphi) // \chi(tanpsi) is the totol deflection angle with parameter tanpsi=-(r √f dphi)/dr
double chi(double cotpsi){ double chi(double tanpsi){
gsl_odeiv2_system sys = {func, jac, 2, NULL}; 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); gsl_odeiv2_driver *d = gsl_odeiv2_driver_alloc_y_new(&sys, gsl_odeiv2_step_rk8pd, 1e-3, 1e-12, 1e-10);
@@ -72,7 +72,7 @@ double chi(double cotpsi){
double y[2] = { double y[2] = {
u0, 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 u0*sqrt(f(R0))/tanpsi
}; };
double phi = 0; double phi = 0;
double phi_max = 8*PI; double phi_max = 8*PI;
@@ -142,14 +142,13 @@ int init(double rmax, double rela_err_limit, int *size, double **x, double **y){
sample.size=1; sample.size=1;
sample.x = malloc(sizeof(double) * 1); sample.x = malloc(sizeof(double) * 1);
sample.y = malloc(sizeof(double) * 1); sample.y = malloc(sizeof(double) * 1);
double cotpsi_min = 1/rmax; sample.x[0] = tanpsi_min;
sample.x[0] = cotpsi_min; sample.y[0] = chi(tanpsi_min);
sample.y[0] = chi(cotpsi_min); double chi_tanpsimax = chi(rmax);
double chi_cotpsimax = chi(cotpsi_max);
refine_interval(cotpsi_min, cotpsi_max, sample.y[0], chi_cotpsimax, rela_err_limit, &sample, 1); refine_interval(tanpsi_min, rmax, sample.y[0], chi_tanpsimax, rela_err_limit, &sample, 1);
sample_push(&sample, cotpsi_max, chi_cotpsimax); sample_push(&sample, rmax, chi_tanpsimax);
//sort //sort

30
src/render.c
View File

@@ -22,15 +22,17 @@ inline double xy_to_b(double x, double y) {
} }
double integrand(double *xy, size_t dim, void *integrand_params_void) { double integrand(double *xy, size_t dim, void *integrand_params_void) {
double cotpsi = 1.0/hypot(xy[0], xy[1]); double tanpsi = hypot(xy[0], xy[1]);
if (cotpsi >= cotpsi_max) return 0; double cotpsi = 1.0/tanpsi;
if (tanpsi <= tanpsi_min) return 0;
Integrand_params *integrand_params = (Integrand_params *)integrand_params_void; 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 theta = gsl_spline_eval(integrand_params->spline_data->spline, tanpsi, integrand_params->spline_data->acc);
double phi = atan2(xy[1], xy[0]); double phi = atan2(xy[1], xy[0]);
double rgb[3] = {0}; double rgb[3] = {0};
double theta_phi[2] = {theta, phi}; double theta_phi[2] = {theta, phi};
integrand_params->system->angle_to_pixel(theta_phi, rgb); integrand_params->system->angle_to_pixel(theta_phi, rgb);
return rgb[integrand_params->c]*fabs(gsl_spline_eval_deriv(integrand_params->spline_data->spline, cotpsi, integrand_params->spline_data->acc))*POW2(sin(theta))*POW2(1+POW2(cotpsi)); //printf("(x,y)=(%g,%g), tanpsi=%g, theta=%g, %g\n", xy[0], xy[1], tanpsi, theta, POW2(sin(theta))*(2+POW2(cotpsi)+POW2(tanpsi)));
return rgb[integrand_params->c]*fabs(gsl_spline_eval_deriv(integrand_params->spline_data->spline, tanpsi, integrand_params->spline_data->acc))*POW2(sin(theta))*(2+POW2(tanpsi)+POW2(cotpsi));
} }
int MC_pixel_render(const System *system, int i, int j, const Spline_data spline_data, double *rgb, int pixel_render_max, double pixel_render_err, gsl_rng *r){ int MC_pixel_render(const System *system, int i, int j, const Spline_data spline_data, double *rgb, int pixel_render_max, double pixel_render_err, gsl_rng *r){
@@ -39,10 +41,13 @@ int MC_pixel_render(const System *system, int i, int j, const Spline_data spline
double y_lu = -j*dw + (system->H)*dw/2; double y_lu = -j*dw + (system->H)*dw/2;
double xl[2] = {x_lu, y_lu - dw}; double xl[2] = {x_lu, y_lu - dw};
double xu[2] = {x_lu + dw, y_lu}; double xu[2] = {x_lu + dw, y_lu};
double err_dw = pixel_render_err*dw*dw;
//printf("integrating for (i,j)=(%d, %d)\n", i, j);
//TODO: fill the MC intergral //TODO: fill the MC intergral
// linear rgb = ∫ rgb(χ, ϕ) sin^2(χ)/sin^2(ψ) |χ'(ψ)| dx dy // linear rgb = ∫ rgb(χ, ϕ) sin^2(χ)/sin^2(ψ) |χ'(ψ)| dx dy
// χ'(ψ) = χ'(cot(ψ)) * cot'(ψ) = - χ'(cot(ψ))/\sin^2(ψ) // χ'(ψ) = χ'(tan(ψ)) * tan'(ψ) = - χ'(tan(ψ))/cos^2(ψ) = - χ'(tan(ψ))(1+tan^2(ψ))
gsl_monte_miser_state *miser_state = gsl_monte_miser_alloc(2); gsl_monte_miser_state *miser_state = gsl_monte_miser_alloc(2);
double color; double color;
@@ -55,17 +60,23 @@ int MC_pixel_render(const System *system, int i, int j, const Spline_data spline
//printf("i=%d, j=%d, rgb=( ", i, j); //printf("i=%d, j=%d, rgb=( ", i, j);
for (int c = 0; c < 3; c++) { for (int c = 0; c < 3; c++) {
integrand_params.c = c; integrand_params.c = c;
if (i==0 && j==0) {
//printf("(x,y) = (%g, %g), intergrand[%d]=%g\n", xl[0], xl[1], c, integrand(xl, 2, &integrand_params));
}
gsl_monte_miser_init(miser_state); gsl_monte_miser_init(miser_state);
gsl_monte_miser_integrate(&F, xl, xu, 2, 100, r, miser_state, &color, &err); gsl_monte_miser_integrate(&F, xl, xu, 2, 100, r, miser_state, &color, &err);
//printf("%g +- %g ", color, err); //printf("%g +- %g ", color, err);
if (err > pixel_render_err) { if (err > err_dw) {
int ncalls = 100*(err/pixel_render_err)*(err/pixel_render_err); //int ncalls = 100*POW2(err/err_dw);
if (ncalls > pixel_render_max) ncalls = pixel_render_max; //printf("(i,j)=(%d,%d), err=%g, err_dw=%g, ncalls=%d\n", i, j, err, err_dw, ncalls);
//if (ncalls < 0 || ncalls > pixel_render_max) ncalls = pixel_render_max;
int ncalls = pixel_render_max;
//printf("(i,j)=(%d,%d), err=%g, err_dw=%g, ncalls=%d\n", i, j, err, err_dw, ncalls);
gsl_monte_miser_integrate(&F, xl, xu, 2, ncalls, r, miser_state, &color, &err); gsl_monte_miser_integrate(&F, xl, xu, 2, ncalls, r, miser_state, &color, &err);
} }
rgb[c] = color/(dw*dw); rgb[c] = color/(dw*dw);
} }
//printf(")\n"); //printf("finished for (i,j)=(%d, %d)\n", i, j);
return 0; return 0;
} }
@@ -108,5 +119,6 @@ int render(System *system, double *buffer, int pixel_render_max, double pixel_re
gsl_interp_accel_free(acc); gsl_interp_accel_free(acc);
gsl_rng_free(r); gsl_rng_free(r);
} }
printf("rendering finished.\n");
return 0; return 0;
} }

14
src/write_png.c
View File

@@ -14,23 +14,11 @@ static inline double linear_to_srgb(double x){
} }
int buffer_normalize_srgb(double *buffer, int W, int H) { int buffer_normalize_srgb(double *buffer, int W, int H) {
double max = 0;
//for (int j = 0; j < H; j++) {
// for (int i = 0; i < W; i++) {
// for (int c = 0; c < 3; c++) {
// if (max < buffer[color_index(i, j, c)]) max = buffer[color_index(i, j, c)];
// }
// }
//}
max = percentile(buffer, W*H*3, cutperc);
printf("max in the buffer: %g\n", max);
if (max == 0) return 1;
#pragma omp parallel for #pragma omp parallel for
for (int j = 0; j < H; j++) { for (int j = 0; j < H; j++) {
for (int i = 0; i < W; i++) { for (int i = 0; i < W; i++) {
for (int c = 0; c < 3; c++) { for (int c = 0; c < 3; c++) {
buffer[color_index(i, j, c)] /= max; buffer[color_index(i, j, c)] = linear_to_srgb(SCALE*buffer[color_index(i, j, c)]);
buffer[color_index(i, j, c)] = linear_to_srgb(buffer[color_index(i, j, c)]);
} }
} }
} }