潜艇航向散布模型分为四种，分别为确定航向模型、概略航向均匀分布模型、概略航向正态分布模型、未知航向均匀分布模型，分别写出c++代码

下面是四种潜艇航向散布模型的C++代码示例：

1. 确定航向模型：

#include <iostream>
#include <cmath>
#include <ctime>
#include <cstdlib>

using namespace std;

const double PI = 3.14159265358979323846;  // 圆周率
const double R = 6371.0;  // 地球半径，单位：千米

double deg2rad(double deg) {
    return deg * PI / 180.0;
}

double rad2deg(double rad) {
    return rad * 180.0 / PI;
}

double distance(double lat1, double lng1, double lat2, double lng2) {
    double radLat1 = deg2rad(lat1);
    double radLat2 = deg2rad(lat2);
    double a = radLat1 - radLat2;
    double b = deg2rad(lng1) - deg2rad(lng2);
    double s = 2 * asin(sqrt(pow(sin(a / 2), 2) + cos(radLat1) * cos(radLat2) * pow(sin(b / 2), 2)));
    s = s * R;
    return s;
}

double bearing(double lat1, double lng1, double lat2, double lng2) {
    double radLat1 = deg2rad(lat1);
    double radLat2 = deg2rad(lat2);
    double radLng1 = deg2rad(lng1);
    double radLng2 = deg2rad(lng2);
    double y = sin(radLng2 - radLng1) * cos(radLat2);
    double x = cos(radLat1) * sin(radLat2) - sin(radLat1) * cos(radLat2) * cos(radLng2 - radLng1);
    double brng = atan2(y, x);
    brng = rad2deg(brng);
    return fmod((brng + 360), 360);
}

int main() {
    srand((unsigned int) time(NULL));

    double lat1 = 35.0;  // 初始纬度
    double lng1 = 135.0;  // 初始经度
    double heading = 0.0;  // 初始航向
    double speed = 10.0;  // 初始速度，单位：千米/小时
    double timeStep = 1.0;  // 时间步长，单位：小时
    double timeLeft = 24.0;  // 总时间，单位：小时
    double lat2, lng2, dist;

    while (timeLeft > 0) {
        // 计算下一个位置和航向
        dist = speed * timeStep;
        lat2 = asin(sin(deg2rad(lat1)) * cos(dist / R) + cos(deg2rad(lat1)) * sin(dist / R) * cos(deg2rad(heading)));
        lng2 = deg2rad(lng1) + atan2(sin(deg2rad(heading)) * sin(dist / R) * cos(deg2rad(lat1)),
                                     cos(dist / R) - sin(deg2rad(lat1)) * sin(lat2));
        lat2 = rad2deg(lat2);
        lng2 = rad2deg(lng2);

        // 输出位置和距离
        dist = distance(lat1, lng1, lat2, lng2);
        cout << "Lat: " << lat2 << ", Lng: " << lng2 << ", Dist: " << dist << endl;

        // 更新位置和时间
        lat1 = lat2;
        lng1 = lng2;
        timeLeft -= timeStep;
    }

    return 0;
}

2. 概略航向均匀分布模型：

#include <iostream>
#include <cmath>
#include <ctime>
#include <cstdlib>

using namespace std;

const double PI = 3.14159265358979323846;  // 圆周率
const double R = 6371.0;  // 地球半径，单位：千米

double deg2rad(double deg) {
    return deg * PI / 180.0;
}

double rad2deg(double rad) {
    return rad * 180.0 / PI;
}

double distance(double lat1, double lng1, double lat2, double lng2) {
    double radLat1 = deg2rad(lat1);
    double radLat2 = deg2rad(lat2);
    double a = radLat1 - radLat2;
    double b = deg2rad(lng1) - deg2rad(lng2);
    double s = 2 * asin(sqrt(pow(sin(a / 2), 2) + cos(radLat1) * cos(radLat2) * pow(sin(b / 2), 2)));
    s = s * R;
    return s;
}

double bearing(double lat1, double lng1, double lat2, double lng2) {
    double radLat1 = deg2rad(lat1);
    double radLat2 = deg2rad(lat2);
    double radLng1 = deg2rad(lng1);
    double radLng2 = deg2rad(lng2);
    double y = sin(radLng2 - radLng1) * cos(radLat2);
    double x = cos(radLat1) * sin(radLat2) - sin(radLat1) * cos(radLat2) * cos(radLng2 - radLng1);
    double brng = atan2(y, x);
    brng = rad2deg(brng);
    return fmod((brng + 360), 360);
}

int main() {
    srand((unsigned int) time(NULL));

    double lat1 = 35.0;  // 初始纬度
    double lng1 = 135.0;  // 初始经度
    double heading = 0.0;  // 初始航向
    double speed = 10.0;  // 初始速度，单位：千米/小时
    double timeStep = 1.0;  // 时间步长，单位：小时
    double timeLeft = 24.0;  // 总时间，单位：小时
    double lat2, lng2, dist;

    while (timeLeft > 0) {
        // 计算随机下一个航向
        heading = fmod(heading + rand() % 180 - 90, 360);

        // 计算下一个位置
        dist = speed * timeStep;
        lat2 = asin(sin(deg2rad(lat1)) * cos(dist / R) + cos(deg2rad(lat1)) * sin(dist / R) * cos(deg2rad(heading)));
        lng2 = deg2rad(lng1) + atan2(sin(deg2rad(heading)) * sin(dist / R) * cos(deg2rad(lat1)),
                                     cos(dist / R) - sin(deg2rad(lat1)) * sin(lat2));
        lat2 = rad2deg(lat2);
        lng2 = rad2deg(lng2);

        // 输出位置和距离
        dist = distance(lat1, lng1, lat2, lng2);
        cout << "Lat: " << lat2 << ", Lng: " << lng2 << ", Dist: " << dist << endl;

        // 更新位置和时间
        lat1 = lat2;
        lng1 = lng2;
        timeLeft -= timeStep;
    }

    return 0;
}

3. 概略航向正态分布模型：

#include <iostream>
#include <cmath>
#include <ctime>
#include <cstdlib>

using namespace std;

const double PI = 3.14159265358979323846;  // 圆周率
const double R = 6371.0;  // 地球半径，单位：千米

double deg2rad(double deg) {
    return deg * PI / 180.0;
}

double rad2deg(double rad) {
    return rad * 180.0 / PI;
}

double distance(double lat1, double lng1, double lat2, double lng2) {
    double radLat1 = deg2rad(lat1);
    double radLat2 = deg2rad(lat2);
    double a = radLat1 - radLat2;
    double b = deg2rad(lng1) - deg2rad(lng2);
    double s = 2 * asin(sqrt(pow(sin(a / 2), 2) + cos(radLat1) * cos(radLat2) * pow(sin(b / 2), 2)));
    s = s * R;
    return s;
}

double bearing(double lat1, double lng1, double lat2, double lng2) {
    double radLat1 = deg2rad(lat1);
    double radLat2 = deg2rad(lat2);
    double radLng1 = deg2rad(lng1);
    double radLng2 = deg2rad(lng2);
    double y = sin(radLng2 - radLng1) * cos(radLat2);
    double x = cos(radLat1) * sin(radLat2) - sin(radLat1) * cos(radLat2) * cos(radLng2 - radLng1);
    double brng = atan2(y, x);
    brng = rad2deg(brng);
    return fmod((brng + 360), 360);
}

double randn(double mu, double sigma) {
    static bool deviateAvailable = false;
    static double storedDeviate;
    double polar, rsquared, var1, var2;

    if (!deviateAvailable) {
        do {
            var1 = 2.0 * ((double) rand() / (double) RAND_MAX) - 1.0;
            var2 = 2.0 * ((double) rand() / (double) RAND_MAX) - 1.0;
            rsquared = var1 * var1 + var2 * var2;
        } while (rsquared >= 1.0 || rsquared == 0.0);

        polar = sqrt(-2.0 * log(rsquared) / rsquared);
        storedDeviate = var1 * polar;
        deviateAvailable = true;
        return mu + sigma * var2 * polar;
    } else {
        deviateAvailable = false;
        return mu + sigma * storedDeviate;
    }
}

int main() {
    srand((unsigned int) time(NULL));

    double lat1 = 35.0;  // 初始纬度
    double lng1 = 135.0;  // 初始经度
    double heading = 0.0;  // 初始航向
    double speed = 10.0;  // 初始速度，单位：千米/小时
    double timeStep = 1.0;  // 时间步长，单位：小时
    double timeLeft = 24.0;  // 总时间，单位：小时
    double lat2, lng2, dist;

    while (timeLeft > 0) {
        // 计算随机下一个航向
        heading = fmod(heading + randn(0.0, 10.0), 360);

        // 计算下一个位置
        dist = speed * timeStep;
        lat2 = asin(sin(deg2rad(lat1)) * cos(dist / R) + cos(deg2rad(lat1)) * sin(dist / R) * cos(deg2rad(heading)));
        lng2 = deg2rad(lng1) + atan2(sin(deg2rad(heading)) * sin(dist / R) * cos(deg2rad(lat1)),
                                     cos(dist / R) - sin(deg2rad(lat1)) * sin(lat2));
        lat2 = rad2deg(lat2);
        lng2 = rad2deg(lng2);

        // 输出位置和距离
        dist =