import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
 
v_sound = 340
 
def calculate_distance(lat1, lon1, lat2, lon2):
    lat_distance = 111.263 * 1000 * abs(lat2 - lat1)
    mid_lat = (lat1 + lat2) / 2
    lon_distance = 97.304 * 1000 * abs(lon2 - lon1) * np.cos(np.radians(mid_lat))
    return np.sqrt(lat_distance ** 2 + lon_distance ** 2)
 
stations = {
    'A': (110.241, 27.204, 824),
    'B': (110.780, 27.456, 727),
    'C': (110.712, 27.785, 742),
    'D': (110.251, 27.825, 850),
    'E': (110.524, 27.617, 786),
    'F': (110.467, 27.921, 678),
    'G': (110.047, 27.121, 575),
}
 
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')
 
for station, data in stations.items():
    lon, lat, elev = data
    ax.scatter([lon], [lat], [elev], color='red', s=100)
 
for station1, data1 in stations.items():
    for station2, data2 in stations.items():
        if station1 < station2:
            lon1, lat1 ,elev1 = data1
            lon2, lat2, elev2 = data2
            distance_m = calculate_distance(lat1, lon1, lat2, lon2)
            ax.plot([lon1, lon2], [lat1, lat2], [elev1, elev2], color='gray', linestyle=':')
        print(station2, data2)
    print(station1, data1)
 
 
 
ax.set_xlabel('Longitude')
ax.set_ylabel('Latitude')
ax.set_zlabel('Elevation (m)')
plt.title('Monitoring Stations and Approximate Distances')
plt.show()

aW1wb3J0IG51bXB5IGFzIG5wCmltcG9ydCBtYXRwbG90bGliLnB5cGxvdCBhcyBwbHQKZnJvbSBtcGxfdG9vbGtpdHMubXBsb3QzZCBpbXBvcnQgQXhlczNECgp2X3NvdW5kID0gMzQwCgpkZWYgY2FsY3VsYXRlX2Rpc3RhbmNlKGxhdDEsIGxvbjEsIGxhdDIsIGxvbjIpOgogICAgbGF0X2Rpc3RhbmNlID0gMTExLjI2MyAqIDEwMDAgKiBhYnMobGF0MiAtIGxhdDEpCiAgICBtaWRfbGF0ID0gKGxhdDEgKyBsYXQyKSAvIDIKICAgIGxvbl9kaXN0YW5jZSA9IDk3LjMwNCAqIDEwMDAgKiBhYnMobG9uMiAtIGxvbjEpICogbnAuY29zKG5wLnJhZGlhbnMobWlkX2xhdCkpCiAgICByZXR1cm4gbnAuc3FydChsYXRfZGlzdGFuY2UgKiogMiArIGxvbl9kaXN0YW5jZSAqKiAyKQoKc3RhdGlvbnMgPSB7CiAgICAnQSc6ICgxMTAuMjQxLCAyNy4yMDQsIDgyNCksCiAgICAnQic6ICgxMTAuNzgwLCAyNy40NTYsIDcyNyksCiAgICAnQyc6ICgxMTAuNzEyLCAyNy43ODUsIDc0MiksCiAgICAnRCc6ICgxMTAuMjUxLCAyNy44MjUsIDg1MCksCiAgICAnRSc6ICgxMTAuNTI0LCAyNy42MTcsIDc4NiksCiAgICAnRic6ICgxMTAuNDY3LCAyNy45MjEsIDY3OCksCiAgICAnRyc6ICgxMTAuMDQ3LCAyNy4xMjEsIDU3NSksCn0KCmZpZyA9IHBsdC5maWd1cmUoZmlnc2l6ZT0oMTAsIDgpKQpheCA9IGZpZy5hZGRfc3VicGxvdCgxMTEsIHByb2plY3Rpb249JzNkJykKCmZvciBzdGF0aW9uLCBkYXRhIGluIHN0YXRpb25zLml0ZW1zKCk6CiAgICBsb24sIGxhdCwgZWxldiA9IGRhdGEKICAgIGF4LnNjYXR0ZXIoW2xvbl0sIFtsYXRdLCBbZWxldl0sIGNvbG9yPSdyZWQnLCBzPTEwMCkKCmZvciBzdGF0aW9uMSwgZGF0YTEgaW4gc3RhdGlvbnMuaXRlbXMoKToKICAgIGZvciBzdGF0aW9uMiwgZGF0YTIgaW4gc3RhdGlvbnMuaXRlbXMoKToKICAgICAgICBpZiBzdGF0aW9uMSA8IHN0YXRpb24yOgogICAgICAgICAgICBsb24xLCBsYXQxICxlbGV2MSA9IGRhdGExCiAgICAgICAgICAgIGxvbjIsIGxhdDIsIGVsZXYyID0gZGF0YTIKICAgICAgICAgICAgZGlzdGFuY2VfbSA9IGNhbGN1bGF0ZV9kaXN0YW5jZShsYXQxLCBsb24xLCBsYXQyLCBsb24yKQogICAgICAgICAgICBheC5wbG90KFtsb24xLCBsb24yXSwgW2xhdDEsIGxhdDJdLCBbZWxldjEsIGVsZXYyXSwgY29sb3I9J2dyYXknLCBsaW5lc3R5bGU9JzonKQogICAgICAgIHByaW50KHN0YXRpb24yLCBkYXRhMikKICAgIHByaW50KHN0YXRpb24xLCBkYXRhMSkKCgoKYXguc2V0X3hsYWJlbCgnTG9uZ2l0dWRlJykKYXguc2V0X3lsYWJlbCgnTGF0aXR1ZGUnKQpheC5zZXRfemxhYmVsKCdFbGV2YXRpb24gKG0pJykKcGx0LnRpdGxlKCdNb25pdG9yaW5nIFN0YXRpb25zIGFuZCBBcHByb3hpbWF0ZSBEaXN0YW5jZXMnKQpwbHQuc2hvdygp