From 8eeb53b9d74754e351eecc5b62c6022b9e4fedaa Mon Sep 17 00:00:00 2001 From: Yingjie Wang Date: Wed, 5 Nov 2025 12:59:14 -0500 Subject: [PATCH] split files --- src/background.c | 94 +----------------------------------------------- src/init.c | 94 ++++++++++++++++++++++++++++++++++++++++++++++++ src/init.h | 1 + 3 files changed, 96 insertions(+), 93 deletions(-) create mode 100644 src/init.c create mode 100644 src/init.h diff --git a/src/background.c b/src/background.c index 28b4795..f2ac737 100644 --- a/src/background.c +++ b/src/background.c @@ -1,100 +1,8 @@ -#include #include -#include -#include #include #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){ diff --git a/src/init.c b/src/init.c new file mode 100644 index 0000000..cfbcd98 --- /dev/null +++ b/src/init.c @@ -0,0 +1,94 @@ +#include +#include +#include +#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); +} diff --git a/src/init.h b/src/init.h new file mode 100644 index 0000000..2a4f0c5 --- /dev/null +++ b/src/init.h @@ -0,0 +1 @@ +double chi(double b);