「蒼き鋼のアルペジオ ‐アルス・ノヴァ‐」におけるC言語/C++ネタについて ref: http://qiita.com/YSRKEN/items/6fc784d4f7f9206537fc

code.cpp

#define _USE_MATH_DEFINES
#include <cmath>
#include <iostream>

/* using宣言 */
using std::cin;
using std::cout;
using std::endl;

/* 定数宣言
 * 「Fundamental Physical Constants --- Complete Listing」
 * http://physics.nist.gov/cuu/Constants/Table/allascii.txt
 */
const double me = 9.1093826e-31;        //電子の質量(kg)
const double epsilon = 8.854187817e-12; //真空の誘電率(F m^-1)
const double h = 6.626070040e-34;       //プランク定数(J s)
const double ele = 1.6021766208e-19;    //電気素量(C)
const double nucle = 6.644657230e-27;   //α線粒子の質量(kg)
const double Z = 2.0;
const double rm =(2 * me * nucle) / (2 * me + nucle) / 2;
const double ac = (ele * ele) / (4 * M_PI * epsilon * rm);
const double MM = 1.0e-14;

/* main関数 */
int main(){
    // 入力
    cout << "r1 between nucleus and electron 1 (MM)？" << endl;
    double r;
    cin >> r;
    cout << "total energy |E| of helium atom (eV) ?" << endl;
    double E;
    cin >> E;
    // 演算
    for(int i = 1; i < 100; i++, r += 1){
        double poten = -(2 * Z * ele * ele) / (4 * M_PI * epsilon * r) + (ele * ele) / (4 * M_PI * epsilon * 2 * r);
        double vya = -E * ele - poten / MM;
        if(vya > 0.0){
            double vyb = sqrt(vya / me);
            double VY = vyb * 1.0e-9;
            double prexx = r;
            double VX = 0.0;
            double WN = 0.0;
            double preyy = 0.0;
            double preVY, preWN, xx;
            do{
                xx = prexx + VX;
                double yy = preyy + VY;
                double vk = VX * VX + VY * VY;
                double leng = sqrt(vk) * MM;
                double wav = h / (rm * sqrt(vk) * 1.0e9);
                double ra = sqrt(prexx * prexx + preyy * preyy) * 1.0e-14;
                double rb = sqrt(4 * prexx * prexx + 2 * preyy * preyy) * 1.0e-14;
                WN += leng / wav;
                preVY = VY;
                preWN = WN;
                prexx *= MM;
                preyy *= MM;
                VX += 1.0e-32 * ac * prexx * (-Z / (ra * ra * ra) + 2.0 / (rb * rb * rb ));
                VY += 1.0e-32 * ac * preyy * (-Z / (ra * ra * ra) + 1.0 / (rb * rb * rb ));
                prexx = xx;
                preyy = yy;
            }while(xx >= 0.0);
            if(fabs(VY) < 0.0001)
                cout << "r1= " << r << " VX= " << VX << " VY= " << VY << " preVY= " << preVY << " midWN= " << (( preWN + WN ) / 2) << endl;
        }
    }
}