commit 0683d721cb72ecfd70553b50138f7202193c808a Author: moxitech Date: Fri Oct 11 01:43:36 2024 +0700 First lib version diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..f46727f --- /dev/null +++ b/.gitignore @@ -0,0 +1,2 @@ +python +cmd \ No newline at end of file diff --git a/go.mod b/go.mod new file mode 100644 index 0000000..bdc90cb --- /dev/null +++ b/go.mod @@ -0,0 +1,5 @@ +module git-enter-code.ru/moxitech/simulation-go + +go 1.22.7 + +require github.com/dhconnelly/rtreego v1.0.0 diff --git a/go.sum b/go.sum new file mode 100644 index 0000000..8fb55e9 --- /dev/null +++ b/go.sum @@ -0,0 +1,2 @@ +github.com/dhconnelly/rtreego v1.0.0 h1:1+V1STGw+zwx7jpvH/fwbeC5w5gZfn+XinARU45oRek= +github.com/dhconnelly/rtreego v1.0.0/go.mod h1:SDozu0Fjy17XH1svEXJgdYq8Tah6Zjfa/4Q33Z80+KM= diff --git a/pkg/simulation.go b/pkg/simulation.go new file mode 100644 index 0000000..13f00a0 --- /dev/null +++ b/pkg/simulation.go @@ -0,0 +1,378 @@ +package simulation + +import ( + "encoding/json" + "fmt" + "math" + "math/rand" + "sync" + + vector "moxitech/drone/simulation-go/pkg/vector" + "moxitech/drone/simulation-go/utils" + + "github.com/dhconnelly/rtreego" // R-Tree для пространственного индекса +) + +type UAV struct { + Path map[int][]float64 + Speed *vector.Vector3 + Freq float64 // Частота передачи + Radius float64 // Радиус покрытия в километрах +} + +type CGS struct { + Position *vector.Vector3 + Radius float64 // Радиус покрытия в километрах + Frequency float64 // Частота передачи в МГц +} + +type Simulation struct { + elevationDefault float64 + distanceDefault float64 + uavs map[int]*UAV + cgs map[int]*CGS + heights map[string]float64 // высоты в формате "lat,lon": высота + idGen func() int + timeline []int + mutex sync.Mutex + rtree *rtreego.Rtree + results []SimulationResult +} + +type SimulationResult struct { + Time int `json:"time"` + UAVs map[int]UAVResult `json:"uavs"` + CGSs map[int]CGSResult `json:"cgs"` +} + +type UAVResult struct { + Position *vector.Vector3 `json:"position"` + Freq float64 `json:"frequency"` + Radius float64 `json:"radius"` + DataRate float64 `json:"data_rate"` +} + +type CGSResult struct { + Position *vector.Vector3 `json:"position"` + Radius float64 `json:"radius"` + Frequency float64 `json:"frequency"` +} + +func NewSimulation() *Simulation { + return &Simulation{ + elevationDefault: 85.0 * (math.Pi / 180), // 15 degrees in radians + distanceDefault: 1000.0, // km + uavs: make(map[int]*UAV), + cgs: make(map[int]*CGS), + heights: make(map[string]float64), + idGen: utils.GetIDGenerator(), + rtree: rtreego.NewTree(2, 25, 50), // Инициализация R-Tree + results: []SimulationResult{}, + } +} + +// Adding CGS (ground station) with coverage radius and frequency +func (s *Simulation) AddCGS(pos []float64, radius float64, frequency float64) { + s.mutex.Lock() + defer s.mutex.Unlock() + cgsID := s.idGen() + cgsVector := new(vector.Vector3).SetOn(0, pos[0], pos[1]) + s.cgs[cgsID] = &CGS{ + Position: cgsVector, + Radius: radius, + Frequency: frequency, + } + s.rtree.Insert(cgsVector) // Добавление в R-Tree +} + +// Adding UAV with its path and coverage radius +func (s *Simulation) AddUAV(path map[int][]float64, speed []float64, radius float64) { + s.mutex.Lock() + defer s.mutex.Unlock() + uavID := s.idGen() + speedVector := &vector.Vector3{X: speed[0], Y: speed[1], Z: speed[2]} + freq := 5030 + rand.Float64()*(5090-5030) // Генерация частоты в диапазоне 5030-5090 МГц + + s.uavs[uavID] = &UAV{ + Path: path, + Speed: speedVector, + Freq: freq, + Radius: radius, + } + + for _, p := range path { + pos := new(vector.Vector3).SetOn(p[0], p[1], p[2]) + s.rtree.Insert(pos) // Добавление позиций в R-Tree + } +} + +// Removing UAV +func (s *Simulation) RemoveUAV(id int) { + s.mutex.Lock() + defer s.mutex.Unlock() + delete(s.uavs, id) +} + +func (s *Simulation) GetPosition(uav *UAV, time int) *vector.Vector3 { + // Determine known time slices + keys := make([]int, 0) + for k := range uav.Path { + keys = append(keys, k) + } + + if _, exists := uav.Path[time]; exists { + return new(vector.Vector3).SetOn(uav.Path[time][0], uav.Path[time][1], uav.Path[time][2]) + } + + // If out of known range, use speed to estimate position + if time < keys[0] || time > keys[len(keys)-1] { + initialTime := keys[len(keys)-1] + initialPosition := &vector.Vector3{ + X: uav.Path[initialTime][1], + Y: uav.Path[initialTime][2], + Z: uav.Path[initialTime][0], + } + deltaTime := float64(time - initialTime) + speedComponent := uav.Speed.Scale(deltaTime) + return initialPosition.Add(speedComponent) + } + + // Otherwise interpolate + keys = append(keys, time) + keys = utils.Sorted(keys) + idx := utils.IndexOf(keys, time) + temp := new(vector.Vector3) + temp.RerpVectors( + &vector.Vector3{X: uav.Path[keys[idx-1]][1], Y: uav.Path[keys[idx-1]][2], Z: uav.Path[keys[idx-1]][0]}, + &vector.Vector3{X: uav.Path[keys[idx+1]][1], Y: uav.Path[keys[idx+1]][2], Z: uav.Path[keys[idx+1]][0]}, + float64(time-keys[idx-1])/float64(keys[idx+1]-keys[idx-1]), + ) + return temp.Copy() +} + +// Setting timeline configuration +func (s *Simulation) SetTimeline(conf []int) { + s.timeline = make([]int, 0) + for i := conf[0]; i < conf[1]; i += conf[2] { + s.timeline = append(s.timeline, i) + } +} + +// Calculating Great Circle Distance (to account for Earth's curvature) +func GreatCircleDistance(v1, v2 *vector.Vector3) float64 { + lat1, lon1 := v1.Y*(math.Pi/180), v1.X*(math.Pi/180) + lat2, lon2 := v2.Y*(math.Pi/180), v2.X*(math.Pi/180) + dlat := lat2 - lat1 + dlon := lon2 - lon1 + a := math.Sin(dlat/2)*math.Sin(dlat/2) + math.Cos(lat1)*math.Cos(lat2)*math.Sin(dlon/2)*math.Sin(dlon/2) + c := 2 * math.Atan2(math.Sqrt(a), math.Sqrt(1-a)) + R := 6371.0 // Радиус Земли в км + return R * c +} + +// FreeSpacePathLoss Calculating Free Space Path Loss (FSPL) +func FreeSpacePathLoss(distance float64, frequency float64) float64 { + return 20*math.Log10(distance*1e3) + 20*math.Log10(frequency) - 20*math.Log10(3e8/frequency) + 32.44 +} + +// Determining Line of Sight (LoS) +// Function to represent Line of Sight (LoS) check considering terrain and Fresnel zones +func (s *Simulation) HasLineOfSight(v1, v2 *vector.Vector3) bool { + // Расчет расстояния и средней точки между v1 и v2 + distance := GreatCircleDistance(v1, v2) + midpoint := &vector.Vector3{ + X: (v1.X + v2.X) / 2, + Y: (v1.Y + v2.Y) / 2, + Z: (v1.Z + v2.Z) / 2, + } + + // Проверка высоты на средней точке для оценки наличия препятствий + key := fmt.Sprintf("%.6f,%.6f", midpoint.X, midpoint.Y) + terrainHeight, exists := s.heights[key] + if exists && midpoint.Z < terrainHeight { + return false + } + + // Проверка первой зоны Френеля + fresnelRadius := FresnelZoneRadius(distance, v1.Length(), v2.Length()) + if midpoint.Z < terrainHeight+fresnelRadius { + return false + } + + return true +} + +// Calculate the radius of the first Fresnel zone +func FresnelZoneRadius(distance, height1, height2 float64) float64 { + wavelength := 3e8 / 5.06e9 // Средняя частота 5060 МГц + return math.Sqrt((wavelength * distance) / (height1 + height2)) +} + +// Calculating connectivity between UAVs and CGSs +func (s *Simulation) CalculateConnectivity() { + for _, uav := range s.uavs { + for time := range s.timeline { + uavPosition := s.GetPosition(uav, time) + if uavPosition == nil { + continue + } + for _, cgs := range s.cgs { + distance := GreatCircleDistance(uavPosition, cgs.Position) + if distance <= uav.Radius && distance <= cgs.Radius && s.HasLineOfSight(uavPosition, cgs.Position) { + fspl := FreeSpacePathLoss(distance, uav.Freq*1e6) // Используем частоту UAV в МГц + snr := 100.0 / (fspl + 1) // Примерный SNR для расчетов + modulation := s.CalculateModulation(distance) + dataRate := s.CalculateDataRate(modulation, snr) + + fmt.Printf("\u0412ремя %d: UAV на позиции (%.2f, %.2f, %.2f) может подключиться к CGS на позиции (%.2f, %.2f, %.2f) с потерями %.2f дБ, модуляция: %s, частота CGS: %.2f МГц, скорость передачи: %.2f Mbps\n", + time, uavPosition.X, uavPosition.Y, uavPosition.Z, + cgs.Position.X, cgs.Position.Y, cgs.Position.Z, fspl, modulation, cgs.Frequency, dataRate) + } + } + s.results = append(s.results, s.CollectResult(time)) + } + } +} + +// Calculate Modulation Scheme Based on Distance +func (s *Simulation) CalculateModulation(distance float64) string { + if distance < 500 { + return "64-QAM" + } else if distance < 1000 { + return "16-QAM" + } else if distance < 2000 { + return "QPSK" + } + return "BPSK" +} +func (s *Simulation) CalculateDataRate(modulation string, snr float64) float64 { + var spectralEfficiency float64 + switch modulation { + case "64-QAM": + spectralEfficiency = 6.0 + case "16-QAM": + spectralEfficiency = 4.0 + case "QPSK": + spectralEfficiency = 2.0 + case "BPSK": + spectralEfficiency = 1.0 + default: + spectralEfficiency = 0.5 + } + + // Подсчет пропускной способности (в Mbps) + bandwidth := 20.0 // Пропускная способность канала в МГц + dataRate := bandwidth * spectralEfficiency * math.Log2(1+snr) + return dataRate +} + +// CollectResult collects the current state of the simulation at a given time +func (s *Simulation) CollectResult(time int) SimulationResult { + uavsResult := make(map[int]UAVResult) + for id, uav := range s.uavs { + position := s.GetPosition(uav, time) + + if position != nil { + uavsResult[id] = UAVResult{ + Position: position, + Freq: uav.Freq, + Radius: uav.Radius, + DataRate: s.CalculateDataRate(s.CalculateModulation(GreatCircleDistance(position, s.cgs[1].Position)), + 100.0/(FreeSpacePathLoss(GreatCircleDistance(position, s.cgs[1].Position), uav.Freq*1e6)+1)), + } + } + + } + + cgsResult := make(map[int]CGSResult) + for id, cgs := range s.cgs { + cgsResult[id] = CGSResult{ + Position: cgs.Position, + Radius: cgs.Radius, + Frequency: cgs.Frequency, + } + } + + return SimulationResult{ + Time: time, + UAVs: uavsResult, + CGSs: cgsResult, + } +} + +// GetResultsJSON returns the results of the simulation as a JSON string +func (s *Simulation) GetResultsJSON() (string, error) { + jsonData, err := json.MarshalIndent(s.results, "", " ") + if err != nil { + return "", err + } + return string(jsonData), nil +} + +// Test function to create example data and run the simulation +func TestSimulation() []byte { + sim := NewSimulation() + + // Add CGSs (Ground Stations) + sim.AddCGS([]float64{55.7558, 37.6176}, 100, 5060) + sim.AddCGS([]float64{56.7558, 37.7176}, 80, 5080) + sim.AddCGS([]float64{34.0522, -118.2437}, 120, 5035) + + // Add UAVs (drones) with paths and speed + sim.AddUAV(map[int][]float64{ + 0: {1000, 55.7558, 37.6176}, + 10: {1100, 56.0, 38.0}, + 20: {2200, 256.5, 78.5}, + }, []float64{10, 0.1, 0.1}, 100) + + sim.AddUAV(map[int][]float64{ + 0: {1050, 55.7558, 37.6176}, + 10: {1100, 56.0, 38.0}, + 20: {1200, 56.5, 38.5}, + }, []float64{100, 0.1, 0.1}, 50) + + sim.AddUAV(map[int][]float64{ + 0: {1010, 55.7558, 37.6176}, + 10: {1100, 56.0, 38.0}, + 20: {1200, 56.5, 38.5}, + }, []float64{10, 0.1, 0.1}, 100) + + sim.AddUAV(map[int][]float64{ + 0: {1500, 40.7128, -74.0060}, + 10: {1600, 41.0, -73.5}, + 20: {1700, 41.5, -73.0}, + }, []float64{12, 0.15, 0.1}, 100) + + sim.AddUAV(map[int][]float64{ + 0: {1400, 34.0522, -118.2437}, + 10: {1500, 34.5, -118.0}, + 20: {1600, 35.0, -117.5}, + }, []float64{15, 0.2, 0.1}, 100) + + sim.AddUAV(map[int][]float64{ + 0: {1300, 48.8566, 2.3522}, + 10: {1350, 49.0, 2.5}, + 20: {1400, 49.5, 3.0}, + }, []float64{8, 0.1, 0.05}, 100) + + // Set timeline configuration [start, end, step] + sim.SetTimeline([]int{0, 50, 1}) + + // Run simulation and print results + fmt.Println("Результаты симуляции:") + for _, t := range sim.timeline { + fmt.Printf("\nВремя %d секунд:\n", t) + for uavID, uav := range sim.uavs { + position := sim.GetPosition(uav, t) + if position != nil { + fmt.Printf("UAV %d на позиции (%.2f, %.2f, %.2f)\n", uavID, position.X, position.Y, position.Z) + } + } + sim.CalculateConnectivity() + } + result, err := json.Marshal(sim.results) + if err != nil { + panic(err) + } + return result +} diff --git a/pkg/vector/vector.go b/pkg/vector/vector.go new file mode 100644 index 0000000..8ce3c7d --- /dev/null +++ b/pkg/vector/vector.go @@ -0,0 +1,67 @@ +package vector + +import ( + "math" + + "github.com/dhconnelly/rtreego" +) + +type Vector3 struct { + X, Y, Z float64 +} + +func (v *Vector3) SetOn(z, x, y float64) *Vector3 { + v.Z = z + v.X = x + v.Y = y + return v +} + +func (v *Vector3) Copy() *Vector3 { + return &Vector3{v.X, v.Y, v.Z} +} + +func (v *Vector3) Subtract(other *Vector3) *Vector3 { + return &Vector3{ + X: v.X - other.X, + Y: v.Y - other.Y, + Z: v.Z - other.Z, + } +} + +func (v *Vector3) Length() float64 { + return math.Sqrt(v.X*v.X + v.Y*v.Y + v.Z*v.Z) +} + +func (v *Vector3) RerpVectors(v1, v2 *Vector3, ratio float64) { + v.X = v1.X + ratio*(v2.X-v1.X) + v.Y = v1.Y + ratio*(v2.Y-v1.Y) + v.Z = v1.Z + ratio*(v2.Z-v1.Z) +} + +func (v *Vector3) AngleTo(other *Vector3) float64 { + cosAlpha := (v.X*other.X + v.Y*other.Y + v.Z*other.Z) / (v.Length() * other.Length()) + return math.Acos(cosAlpha) +} + +// Implementing rtreego.Spatial interface for Vector3 +func (v *Vector3) Bounds() *rtreego.Rect { + point := rtreego.Point{v.X, v.Y} + rect, _ := rtreego.NewRect(point, []float64{0.001, 0.001}) + return rect +} +func (v *Vector3) Scale(factor float64) *Vector3 { + return &Vector3{ + X: v.X * factor, + Y: v.Y * factor, + Z: v.Z * factor, + } +} + +func (v *Vector3) Add(other *Vector3) *Vector3 { + return &Vector3{ + X: v.X + other.X, + Y: v.Y + other.Y, + Z: v.Z + other.Z, + } +} diff --git a/utils/helpers.go b/utils/helpers.go new file mode 100644 index 0000000..44292bb --- /dev/null +++ b/utils/helpers.go @@ -0,0 +1,37 @@ +package utils + +import "sort" + +// Автоинкрементер ID +func GetIDGenerator() func() int { + id := 0 + return func() int { + id++ + return id + } +} + +// Helper function сортирует срез массивов +func Sorted(values []int) []int { + sort.Ints(values) + return values +} + +// Helper function ищет индекс по значению в слайсе +func IndexOf(values []int, target int) int { + for i, v := range values { + if v == target { + return i + } + } + return -1 +} + +func FindIndex(slice []int, predicate func(int) bool) int { + for i, v := range slice { + if predicate(v) { + return i + } + } + return -1 +}