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

fix: fix intergrad

Browse Source
This commit is contained in:
Yingjie Wang
2025-12-07 22:41:59 -05:00
parent bbbebf1fdf
commit 480c49dcf0
2 changed files with 11 additions and 5 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
src/common.h
View File

@@ -5,6 +5,8 @@
#include <gsl/gsl_errno.h>
#include <math.h>
#define POW2(x) (x*x)
#define FLAT 0
#define PI 3.1415926535897932384626433832795028841971693993751058
#define Rs 1

14
src/render.c
View File

@@ -6,6 +6,7 @@
#include <math.h>
#include <stdio.h>
#include "render.h"
#include "common.h"
typedef struct {
const System *system;
@@ -28,7 +29,8 @@ double integrand_r(double *xy, size_t dim, void *integrand_params_void) {
double rgb[3] = {0};
double theta_phi[2] = {theta, phi};
integrand_params->system->angle_to_pixel(theta_phi, rgb);
return rgb[0]*fabs(gsl_spline_eval_deriv(integrand_params->spline_data->spline, cotpsi, integrand_params->spline_data->acc));
return rgb[0]*fabs(gsl_spline_eval_deriv(integrand_params->spline_data->spline, cotpsi, integrand_params->spline_data->acc))*POW2(sin(theta))*POW2(1+POW2(cotpsi));
return rgb[0];
}
double integrand_g(double *xy, size_t dim, void *integrand_params_void) {
@@ -40,7 +42,8 @@ double integrand_g(double *xy, size_t dim, void *integrand_params_void) {
double rgb[3] = {0};
double theta_phi[2] = {theta, phi};
integrand_params->system->angle_to_pixel(theta_phi, rgb);
return rgb[1]*fabs(gsl_spline_eval_deriv(integrand_params->spline_data->spline, cotpsi, integrand_params->spline_data->acc));
return rgb[1]*fabs(gsl_spline_eval_deriv(integrand_params->spline_data->spline, cotpsi, integrand_params->spline_data->acc))*POW2(sin(theta))*POW2(1+POW2(cotpsi));
return rgb[1];
}
double integrand_b(double *xy, size_t dim, void *integrand_params_void) {
@@ -52,7 +55,8 @@ double integrand_b(double *xy, size_t dim, void *integrand_params_void) {
double rgb[3] = {0};
double theta_phi[2] = {theta, phi};
integrand_params->system->angle_to_pixel(theta_phi, rgb);
return rgb[2]*fabs(gsl_spline_eval_deriv(integrand_params->spline_data->spline, cotpsi, integrand_params->spline_data->acc));
return rgb[2]*fabs(gsl_spline_eval_deriv(integrand_params->spline_data->spline, cotpsi, integrand_params->spline_data->acc))*POW2(sin(theta))*POW2(1+POW2(cotpsi));
return rgb[2];
}
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){
@@ -63,7 +67,7 @@ int MC_pixel_render(const System *system, int i, int j, const Spline_data spline
double xu[2] = {x_lu + dw, y_lu};
//TODO: fill the MC intergral
// linear rgb = ∫ rgb(χ, ϕ) sin^2(χ) |χ'(ψ)| dx dy
// linear rgb = ∫ rgb(χ, ϕ) sin^2(χ)/sin^2(ψ) |χ'(ψ)| dx dy
// χ'(ψ) = χ'(cot(ψ)) * cot'(ψ) = - χ'(cot(ψ))/\sin^2(ψ)
gsl_monte_miser_state *miser_state = gsl_monte_miser_alloc(2);
@@ -94,7 +98,7 @@ int MC_pixel_render(const System *system, int i, int j, const Spline_data spline
if (ncalls > pixel_render_max) ncalls = pixel_render_max;
gsl_monte_miser_integrate(&F, xl, xu, 2, ncalls, r, miser_state, &color, &err);
}
rgb[c] = color;
rgb[c] = color/(dw*dw);
}
//printf(")\n");