「蒼き鋼のアルペジオ ‐アルス・ノヴァ‐」におけるC言語/C++ネタについて

#概要
　レンタル屋で蒼き鋼のアルペジオ ‐アルス・ノヴァ‐を借りて観ました。とても面白かったのですが、霧の艦隊キャラの周りにぐるぐる回っているのはどう見てもCかC++ですよねこれ。
　最初は全部手打ちしようかと思いましたが、幸いこれらのソースがあるサイトからのコピペ(許可は取っていると思うが)だそうなので、リンクを貼るだけで纏められて実に楽ちんです。
　ということで、アニメ版アルペにおけるソースコードをざざっと紹介。

#各話リスト
##第1話 航路を持つ者
　冒頭でイオナらがナガラを沈める際にイオナの周りを回っていたコードがこちら。
　　サンプル C 言語コンピュータープログラム（ヘリウム原子軌道計算）
　　Sample list 6-4　任意の型の座標を格納できるPointクラス
　この2つは、イオナが千早群像を連れて横須賀を脱出する際にも流れていました。

##第2話 嵐の中へ
　冒頭でタカオの周囲を流れていたのは……このコードですね。
　　Sample list 6-4　任意の型の座標を格納できるPointクラス
　で、潜水艦内での会話シーンの後、再びタカオの画面で他に流れていたのはこれでした。
　　サンプル C 言語コンピュータープログラム（ヘリウム原子軌道計算）
　……あれ、つまり使い回し？　この後のイオナのシーンでも同じだったし……。

##第3話以降
　第12話まで一通り確認しましたが、上記の2コード以外は見受けられませんでした。ただ、見逃しがないとは言い切れませんので、気づいた方はコメントまでお願いします。
　どうせならLinuxのコードとか色々使ったらよかったんじゃないかと思いましたが、アニメにGPLのコード使うのってセーフなのかがよく分かりません……。

##劇場版 蒼き鋼のアルペジオ -アルス・ノヴァ- DC
　総集編＋αということで、新規コードがあるんじゃないかと期待していました。
　そうだよ期待してたんだよ！なんでヒエイも上2つのコードしか無いんだよ！！
　もう一つの劇場版は……円盤出てから検証を書くかな。

#おまけ
　あくまで個人的な趣味なのですが、C言語の方は凄く書き直したくなりません？
　というわけで、C++風に書き直してみました。
(注意：サンプルデータとしては、r = 3060.0, E = 79.0を使いましょう。
詳しい意味は、当該サイトをご覧ください)

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;
		}
	}
}