Txgy.EWS.Client/Txgy.EWS.Client.FocalMechanism/Core/PolygonFactory.cs

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using Txgy.EWS.Client.FocalMechanism.Model;
using Txgy.Microseismic.BaseLib.Models;
using FmGrid = Txgy.EWS.Client.FocalMechanism.Model.FmGrid;
using FMStation = Txgy.EWS.Client.FocalMechanism.Model.FMStation;

namespace Txgy.EWS.Client.FocalMechanism.Core
{
    /// <summary>
    /// 张性
    /// </summary>
    public class PolygonFactory : Factory
    {
        public double[] posAve { get; set; }
        public double[] negAve { get; set; }

        public PolygonFactory(FmGrid g) : base(g)
        {
            posAve = new double[180];
            negAve = new double[180];
        }

        public void initAve(FMMap map)
        {
            for (int i = 0; i < 180; i++)
            {
                int posC = 0;
                for (int j = 0; j < map.posStation.Count; j++)
                {
                    int fj = grid.StationList.ToList<FMStation>().FindIndex(fi => fi.Name == map.posStation[j]);
                    if (fj == -1) continue;
                    //if (grid.stations[map.posStation[j]] == null) continue;
                    //posAve[i] += grid.distance[i, map.posStation[j]];
                    posAve[i] += grid.distance[i, fj];
                    posC++;
                }
                posAve[i] /= posC;

                int negC = 0;
                for (int j = 0; j < map.negStation.Count; j++)
                {
                    int fj = grid.StationList.ToList<FMStation>().FindIndex(fi => fi.Name == map.negStation[j]);
                    if (fj == -1) continue;
                    negAve[i] += grid.distance[i, fj];
                    //if (grid.stations[map.negStation[j]] == null) continue;
                    //negAve[i] += grid.distance[i, map.negStation[j]];
                    negC++;
                }
                negAve[i] /= negC;
            }
        }

        /// <summary>
        /// 判断是否为张性
        /// </summary>
        /// <param name="map"></param>
        /// <returns></returns>
        public ComputationResult generateResult(FMMap map)
        {
            ComputationResult r = null;
            Polygon p = null;
            if (map.posStation.Count == 1)
            {
                p = onePointCircle(map, map.posStation[0], map.negStation);
            }
            else if(map.negStation.Count == 1){
                p = onePointCircle(map, map.negStation[0], map.posStation);
            }
            else p = this.computePolygon(map);
            if (p != null)
            {
                r = new ComputationResult();
                r.graph = p;
                r.Direction = p.axis.angle < 180 ? p.axis.angle : p.axis.angle - 180;
                r.FocalType = FocalMechanismType.tensile;
                map.T = FocalMechanismType.tensile;
            }
            return r;
        }

        private Polygon onePointCircle(FMMap map, string point, List<string> stations)
        {
            double diff = 0;
            int angle = -1;
            int pi = grid.StationList.ToList<FMStation>().FindIndex(fi => fi.Name == point);
            for (int i = 0; i < 180; i++)
            {
                double ceiling = int.MaxValue;
                double floor = int.MinValue;
                
                if (pi == -1) continue;
                for (int j = 0; j < stations.Count; j++)
                {
                    int fj = grid.StationList.ToList<FMStation>().FindIndex(fi => fi.Name == stations[j]);
                    if (fj == -1) continue;
                    //if (grid.stations[stations[j]] == null) continue;
                    if (grid.distance[i, fj] > grid.distance[i, pi])
                    {
                        ceiling = Math.Min(ceiling, grid.distance[i, fj]);
                    }
                    else
                    {
                        floor = Math.Max(floor, grid.distance[i, fj]);
                    }
                }

                if (ceiling - floor > diff)
                {
                    angle = i;
                    diff = ceiling - floor;
                }
            }
            Polygon p = new Polygon();
            p.addPoint(this.grid.StationList[pi]);
            p.axis = new Line(angle, grid.distance[angle, pi]);
            return p;
        }

        private Polygon computePolygon(FMMap map)
        {
            this.initAve(map);
            Polygon posp = new Polygon();
            Polygon negp = new Polygon();
            bool posCorrect = true;
            bool negCorrect = true;
            for (int i = 0; i < map.posStation.Count; i++)
            {
                int fj = grid.StationList.ToList<FMStation>().FindIndex(fi => fi.Name == map.posStation[i]);
                if (fj == -1) continue;
                posp.addPoint(this.grid.StationList[fj]);
            }

            for (int i = 0; i < map.negStation.Count; i++)
            {
                int fj = grid.StationList.ToList<FMStation>().FindIndex(fi => fi.Name == map.negStation[i]);
                if (fj == -1) continue;
                if (posp.contains(this.grid.StationList[fj]))
                {
                    posCorrect = false;
                    break;
                }
            }
            if (posCorrect)
            {
                posp.axis = this.computePosAxis(map);
                return posp;
            }

            for (int i = 0; i < map.negStation.Count; i++)
            {
                int fj = grid.StationList.ToList<FMStation>().FindIndex(fi => fi.Name == map.negStation[i]);
                if (fj == -1) continue;
                negp.addPoint(this.grid.StationList[fj]);
            }

            for (int i = 0; i < map.posStation.Count; i++)
            {
                int fj = grid.StationList.ToList<FMStation>().FindIndex(fi => fi.Name == map.posStation[i]);
                if (fj == -1) continue;
                if (negp.contains(this.grid.StationList[fj]))
                {
                    negCorrect = false;
                    break;
                }
            }
            if (negCorrect)
            {
                negp.axis = this.computeNegAxis(map);
                return negp;
            }
            return null;
        }

        private Line computePosAxis(FMMap map)
        {
            int angle = 0;
            double Dis = 0;
            double minDis = int.MaxValue;
            for (int i = 0; i < 180; i++)
            {
                double disSum = 0;
                for (int j = 0; j < map.posStation.Count; j++)
                {
                    int fj = grid.StationList.ToList<FMStation>().FindIndex(fi => fi.Name == map.posStation[j]);
                    if (fj == -1) continue;
                    disSum += Math.Abs(posAve[i] - this.grid.distance[i, fj]);
                }
                if (disSum < minDis)
                {
                    minDis = disSum;
                    angle = i;
                    Dis = posAve[i];
                }
            }
            return new Line(angle, Dis);
        }

        private Line computeNegAxis(FMMap map)
        {
            int angle = 0;
            double Dis = 0;
            double minDis = int.MaxValue;
            for (int i = 0; i < 180; i++)
            {
                double disSum = 0;
                for (int j = 0; j < map.negStation.Count; j++)
                {
                    int fj = grid.StationList.ToList<FMStation>().FindIndex(fi => fi.Name == map.negStation[j]);
                    if (fj == -1) continue;
                    //if (this.grid.stations[map.negStation[j]] == null) continue;
                    disSum += Math.Abs(negAve[i] - this.grid.distance[i, fj]);
                }
                if (disSum < minDis)
                {
                    minDis = disSum;
                    angle = i;
                    Dis = negAve[i];
                }
            }
            return new Line(angle, Dis);
        }
    }
}