碼農的日誌

一

intersection裡面會有負數開平方的問題，因為兩個點離得太近了。

一個trivial的解決辦法：

算r3的時候：

pp = cos(z_2 - z_1);

if (std::fabs(pp - 1.00) < 1e-5) r3 = std::fabs(r1 - r2);

else r3 = sqrt(r1*r1 + r2*r2 - 2 * r1*r2*pp);

算xd_3的時候：

if (std::fabs(pp - 1.00) < 1e-5) x_d[3] = 0.5*z_1 + 0.5*z_2;

else ...

二

博後師兄寫raytrace裡面定義的一堆常數有個命名成了d1正好和defpar重了，改成constd1了，還好後面也只用了一次，沒弄混。

三

raytrace最後一步算cosem的時候用新的代碼正好和原來差了個負號

redshift_KRZ(spin,d1,xem[1], const1, gfactor);

// /*Non Kerr PRD 90, 064002 (2014) Eq. 34*/

// cosem = carter*gfactor/sqrt(xem[1]*xem[1]);

long double g[4][4];

metric_KRZ(spin, d1, xem[1], Pi / 2.0, g);

cosem = -gfactor*sqrt(g[2][2])*kth;

四

VS里用不了庫函數的max於是隨手寫了也不知道會不會慢，懷念起了高中時光

long double mymax(long double a, long double b) {

if (a > b) return a;

else return b;

}

五

看著上角那蜜汁凸起，舊代碼不用數值算metric的時候也有一模一樣的凸起。。。

要算精確的話又要算這一坨東西的求導了。。。

六

visual studio的麻煩真的好多。。。

sprintf(filename_o,"photon/photons4trf_a%.05Le.i%.02Le.e_%.02Le.a13_%.02Le.a22_%.02Le.a52_%.02Le.dat",spin,iobs_deg,epsi3,a13,a22,a52);

注意在visual studio這裡改成

sprintf(filename_o,"photons4trf_a%.05Le.i%.02Le.e_%.02Le.a13_%.02Le.a22_%.02Le.a52_%.02Le.dat",spin,iobs_deg,epsi3,a13,a22,a52);

七

2018-4-8發現有一個求導少乘了個2。。。。

八

2018-4-12

新代碼還是有和以前一樣的問題，raytrace會在度規不可逆的地方死循環，加這麼幾句吧。。

long double gg;

long double myg[4][4];

metric(r, th, myg);

gg = myg[0][0] * myg[3][3] - myg[0][3] * myg[3][0];

if (std::fabs(gg) < 1e-9) stop_integration = 4;

九

2018-4-13

大概case<20的時候會波動很大。。。不過估計沒什麼事吧，畢竟半徑那麼大的地方肯定已經幾乎沒反射譜的貢獻了吧

推薦閱讀：