C,
pasted
on Feb 19:
|
|
#define _USE_MATH_DEFINES
#include <stdio.h>
#include <math.h>
/* 無次元化レナード・ジョーンズ・ポテンシャル */
double ndLJPotential(double x){
/* y = x^(-6) */
double y = 1 / x; y *= y; y *= y*y;
/* 4(x^(-12) - x^(-6) */
return 4 * ( y - 1 ) * y;
}
/* 転回点を求める */
double getXin(double energy, double precision){
double x_min = 1.12; /* x_min = r_min / a ≒ 1.12 */
double dx = x_min / 10;
double x1 = x_min;
double x2 = x_min - dx;
/* 符号の変わり目を挟むまで移動 */
while(ndLJPotential(x2) - energy < 0){
x1 = x2;
x2 -= dx;
}
/* 精度が良くなるまで繰り返す */
while(x1-x2 > precision){
double x = (x1+x2)/2;
if(ndLJPotential(x) - energy > 0){
x2 = x;
}else{
x1 = x;
}
}
return (x1+x2) / 2;
}
/* 転回点を求める */
double getXout(double energy, double precision){
double x_min = 1.12; /* x_min = r_min / a ≒ 1.12 */
double dx = x_min / 10;
double x1 = x_min;
double x2 = x_min + dx;
/* 符号の変わり目を挟むまで移動 */
while(ndLJPotential(x2) - energy < 0){
x1 = x2;
x2 += dx;
}
/* 精度が良くなるまで繰り返す */
while(x2-x1 > precision){
double x = (x1+x2)/2;
if(ndLJPotential(x) - energy > 0){
x2 = x;
}else{
x1 = x;
}
}
return (x1+x2) / 2;
}
/* 量子化されたエネルギーを求める */
double getQuantizedEnergy(double energy, double gamma, int N, double precision){
int i;
double xin = getXin(energy, precision);
double xout = getXout(energy, precision);
/* 合成シンプソン1/3則で積分 */
/* 最初と最後は積分範囲の決め方から0になるので無視 */
double sum = 0.0;
double dx = (xout - xin) / N;
/*sum += sqrt(energy - ndLJPotential(xin));*/
for(i = 1; i < N-1; i++){
if(i%2 == 1){
sum += 4 * sqrt(energy - ndLJPotential(xin + i*dx));
}else{
sum += 2 * sqrt(energy - ndLJPotential(xin + i*dx));
}
}
/*sum += sqrt(energy - ndLJPotential(xout));*/
sum *= dx / 3.0;
return 2 * gamma * sum;
}
/* 量子数を求める */
double getQuantumNumber(double energy, double gamma, int N, double precision){
double s = getQuantizedEnergy(energy, gamma, N, precision);
return s / (2*M_PI) - 0.5;
}
/* エネルギー準位を求める */
double getEnergyLevel(int n, double e1, double e2, double gamma, int N, double x_precision, double energy_precision, int* pitr){
int itr = 0;
double f1 = getQuantizedEnergy(e1, gamma, N, x_precision) - 2*M_PI*(n+0.5);
double f2 = getQuantizedEnergy(e2, gamma, N, x_precision) - 2*M_PI*(n+0.5);
while(fabs(f2) > energy_precision){
double e = (e1*f2-e2*f1)/(f2-f1);
e1 = e2;
e2 = e;
f1 = f2;
f2 = getQuantizedEnergy(e2, gamma, N, x_precision) - 2*M_PI*(n+0.5);
itr++;
}
*pitr = itr;
return e2;
}
/* 全てのエネルギー順位を求める */
void printEnergyLevels(double max_energy, double init_energy, double gamma, int N, double x_precision, double energy_step, double energy_precision){
int i;
int itr;
double xin, xout;
double energy = init_energy + energy_step;
int n = (int)(getQuantumNumber(max_energy, gamma, N, x_precision));
for(i = 1; i <= n; i++){
energy = getEnergyLevel(i, energy, energy+energy_step, gamma, N, x_precision, energy_precision, &itr);
xin = getXin(energy, x_precision);
xout = getXout(energy, x_precision);
printf("s%2d=%f, x_in=%f, x_out=%f, iteration=%d\n", i, energy, xin, xout, itr);
}
}
int main(){
/* step1 */
/* printf("xin = %f, xout = %f\n", getXin(-0.75, 10e-8), getXout(-0.75, 10e-8)); */
/* step2 */
/* printf("int = %f\n", getQuantizedEnergy(-0.75, 100, 100, 10e-8)); */
/* step2 */
/* printf("qnb = %f\n", getQuantumNumber(-0.0005, 100, 100, 10e-8)); */
printEnergyLevels(
-0.0005, /* 最大エネルギー */
-1, /* 初期エネルギー */
100, /* γ */
100, /* 積分の分割数 */
10e-8, /* x_in, x_outの精度 */
0.001, /* エネルギーの刻み幅 */
0.0005) /* エネルギーの精度 */
;
return 0;
}
|
Output:
|
|
s 1=-0.847207, x_in=1.062404, x_out=1.219145, iteration=3
s 2=-0.754142, x_in=1.049601, x_out=1.258183, iteration=3
s 3=-0.667862, x_in=1.040474, x_out=1.295183, iteration=3
s 4=-0.588154, x_in=1.033444, x_out=1.331909, iteration=3
s 5=-0.514801, x_in=1.027803, x_out=1.369284, iteration=3
s 6=-0.447582, x_in=1.023164, x_out=1.407946, iteration=3
s 7=-0.386273, x_in=1.019287, x_out=1.448431, iteration=3
s 8=-0.330643, x_in=1.016016, x_out=1.491248, iteration=3
s 9=-0.280456, x_in=1.013237, x_out=1.536932, iteration=3
s10=-0.235470, x_in=1.010871, x_out=1.586080, iteration=3
s11=-0.195438, x_in=1.008854, x_out=1.639383, iteration=3
s12=-0.160104, x_in=1.007137, x_out=1.697674, iteration=3
s13=-0.129208, x_in=1.005681, x_out=1.761981, iteration=3
s14=-0.102482, x_in=1.004454, x_out=1.833602, iteration=3
s15=-0.079651, x_in=1.003428, x_out=1.914219, iteration=3
s16=-0.060432, x_in=1.002580, x_out=2.006073, iteration=4
s17=-0.044537, x_in=1.001889, x_out=2.112208, iteration=4
s18=-0.031670, x_in=1.001336, x_out=2.236944, iteration=4
s19=-0.021527, x_in=1.000905, x_out=2.386632, iteration=4
s20=-0.013799, x_in=1.000578, x_out=2.571101, iteration=4
s21=-0.008168, x_in=1.000341, x_out=2.806598, iteration=4
s22=-0.004311, x_in=1.000180, x_out=3.122519, iteration=4
s23=-0.001900, x_in=1.000079, x_out=3.579623, iteration=4
s24=-0.000603, x_in=1.000025, x_out=4.334576, iteration=4
|
|