用python 写出程序实现：构建坐标系图，用户输入A点（A点为地面监测站）经纬度，海拔，图上显示在相应的位置并用五角星表示，然后输入B点经纬度，海拔B点在图上相应的位置显示并用三角号表示，A点的检测距离是200km，方位辐射角度是40度，俯仰角度是15度，并在图上用扇形阴影面表示，

好的，这是一个比较复杂的程序，需要用到许多库和数学公式，我将逐步给出程序实现的步骤和代码示例。

首先，我们需要导入所需的库，包括math、numpy、matplotlib等。其中，math库提供了一些常用的数学函数，numpy库提供了高性能的数组运算，matplotlib库用于绘制图形。

import math
import numpy as np
import matplotlib.pyplot as plt

接下来，我们需要定义一些常量和参数，包括地球半径、经纬度转换参数、A点和B点的经纬度和海拔等信息。

# 地球半径，单位为千米
EARTH_RADIUS = 6371.0

# 经纬度转换参数
DEG2RAD = math.pi / 180.0
RAD2DEG = 180.0 / math.pi

# A点经纬度和海拔
ALAT, ALON, AALT = 30.0, 120.0, 50.0

# B点经纬度和海拔
BLAT, BLON, BALT = 31.0, 121.0, 100.0

# A点检测距离，单位为千米
DISTANCE = 200.0

# 方位辐射角度，单位为度
AZIMUTH_ANGLE = 40.0

# 俯仰角度，单位为度
ELEVATION_ANGLE = 15.0

接下来，我们需要定义一些辅助函数，包括经纬度转换函数、距离计算函数、方位角计算函数、俯仰角计算函数等。

def haversine(lat1, lon1, lat2, lon2):
    """
    计算两点之间的距离
    :param lat1: 第一个点的纬度，单位为度
    :param lon1: 第一个点的经度，单位为度
    :param lat2: 第二个点的纬度，单位为度
    :param lon2: 第二个点的经度，单位为度
    :return: 两点之间的距离，单位为千米
    """
    dlat = math.radians(lat2 - lat1)
    dlon = math.radians(lon2 - lon1)
    a = math.sin(dlat / 2) ** 2 + math.cos(math.radians(lat1)) * math.cos(math.radians(lat2)) * math.sin(dlon / 2) ** 2
    c = 2 * math.atan2(math.sqrt(a), math.sqrt(1 - a))
    return EARTH_RADIUS * c

def bearing(lat1, lon1, lat2, lon2):
    """
    计算两点之间的方位角
    :param lat1: 第一个点的纬度，单位为度
    :param lon1: 第一个点的经度，单位为度
    :param lat2: 第二个点的纬度，单位为度
    :param lon2: 第二个点的经度，单位为度
    :return: 两点之间的方位角，单位为度
    """
    dlat = math.radians(lat2 - lat1)
    dlon = math.radians(lon2 - lon1)
    y = math.sin(dlon) * math.cos(math.radians(lat2))
    x = math.cos(math.radians(lat1)) * math.sin(math.radians(lat2)) - math.sin(math.radians(lat1)) * \
        math.cos(math.radians(lat2)) * math.cos(dlon)
    return math.atan2(y, x) * RAD2DEG

def elevation(lat1, lon1, alt1, lat2, lon2, alt2):
    """
    计算两点之间的俯仰角
    :param lat1: 第一个点的纬度，单位为度
    :param lon1: 第一个点的经度，单位为度
    :param alt1: 第一个点的海拔，单位为千米
    :param lat2: 第二个点的纬度，单位为度
    :param lon2: 第二个点的经度，单位为度
    :param alt2: 第二个点的海拔，单位为千米
    :return: 两点之间的俯仰角，单位为度
    """
    dist = haversine(lat1, lon1, lat2, lon2)
    elv = math.atan2(alt2 - alt1, dist)
    return elv * RAD2DEG

接下来，我们需要根据A点的经纬度和海拔，计算出A点在图上的坐标，并用五角星表示。

# 计算A点在图上的坐标
ax, ay = np.array([ALON, ALAT]) * DEG2RAD
ax, ay = EARTH_RADIUS * math.cos(ay) * math.sin(ax), EARTH_RADIUS * math.sin(ay)
ax, ay = ax / 1000.0, ay / 1000.0

# 绘制A点
plt.plot(ax, ay, marker='*', markersize=10, markerfacecolor='yellow', markeredgecolor='black')

接下来，我们需要根据B点的经纬度和海拔，计算出B点在图上的坐标，并用三角号表示。

# 计算B点在图上的坐标
bx, by = np.array([BLON, BLAT]) * DEG2RAD
bx, by = EARTH_RADIUS * math.cos(by) * math.sin(bx), EARTH_RADIUS * math.sin(by)
bx, by = bx / 1000.0, by / 1000.0

# 绘制B点
plt.plot(bx, by, marker='^', markersize=10, markerfacecolor='green', markeredgecolor='black')

接下来，我们需要根据A点的检测距离、方位辐射角度和俯仰角度，计算出扇形阴影面的边界，并在图上绘制出来。

# 计算扇形阴影面的边界
theta = np.linspace(-AZIMUTH_ANGLE / 2, AZIMUTH_ANGLE / 2, 50) * DEG2RAD
phi = np.linspace(0, ELEVATION_ANGLE, 50) * DEG2RAD
r = DISTANCE * 1000.0

x = r * np.outer(np.sin(theta), np.cos(phi)) + ax
y = r * np.outer(np.cos(theta), np.cos(phi)) + ay
z = r * np.outer(np.ones(np.size(theta)), np.sin(phi))

# 绘制扇形阴影面
ax = plt.gca(projection='3d')
ax.plot_surface(x, y, z, color='gray', alpha=0.5)

最后，我们需要设置一些图形参数，如坐标轴范围、标题、标签等，并显示图形。

# 设置坐标轴范围
plt.xlim([min(ax, bx) - 1, max(ax, bx) + 1])
plt.ylim([min(ay, by) - 1, max(ay, by) + 1])

# 设置标题和标签
plt.title('Coordinate System')
plt.xlabel('Longitude (km)')
plt.ylabel('Latitude (km)')

# 显示图形
plt.show()

完整代码如下：