WindStressPRM.cs

using System;
using System.Collections.Generic;

namespace WindStressPRM
{
    /// <summary>
    /// Index class for raster lists in list
    /// Индекс-класс для растровых массивов
    /// </summary>
    public class Index
    {
        /// <summary>
        /// Outer index
        /// Внешний индекс
        /// </summary>
        public int Row { get; set; }
        /// <summary>
        /// Inner index
        /// Внутренний индекс
        /// </summary>
        public int Col { get; set; }
        /// <summary>
        /// designated constructor
        /// </summary>
        /// <param name="Row">row index</param>
        /// <param name="Col">column index</param>
        public Index(int Row, int Col)
        {
            if (Row <= -1 || Col <= -1) {
                throw new System.ArgumentOutOfRangeException("Index must be initialized with nonegative integer value");
            }
            this.Row = Row;
            this.Col = Col;
        }
    }
    /// <summary>
    /// Coordinate pair 
    /// </summary>
    public class Coordinate
    {
        /// <summary>
        /// easting coordinate
        /// Координата по широте
        /// </summary>
        public double X { get; private set;}
        /// <summary>
        /// northing coordinate
        /// Координата по долготе
        /// </summary>
        public double Y { get; private set; }
        /// <summary>
        /// designated constructor
        /// Основной конструктор
        /// </summary>
        /// <param name="X"></param>
        /// <param name="Y"></param>
        public Coordinate(double X, double Y)
        {
            this.X = X;
            this.Y = Y;
            if (!(this.CheckValue()))
            {
                throw new System.ArgumentException("Passed coordinates are not valid!");
            }
        }
        /// <summary>
        /// Проверка на валидность значений
        /// </summary>
        /// <returns></returns>
        public bool CheckValue()
        {
            return !Double.IsNaN(X) && !Double.IsInfinity(X) && !Double.IsNaN(Y) && !Double.IsInfinity(Y);
        }
    }
    /// <summary>
    /// Cell obj for regular prognostic wind field
    /// Объект ячейки для поля прогностического ветра на высоте 10м на регулярной сетке
    /// </summary>
    public class PrognosticCell
    {
        /// <summary>
        /// U - component of wind velocity, m/s
        /// U - компонента скорости ветра, м/с
        /// </summary>
        public double VelocityX { get; private set; }
        /// <summary>
        /// V - component of wind velocity, m/s
        /// V - компонента скорости ветра, м/с 
        /// </summary>
        public double VelocityY { get; private set; }
        /// <summary>
        /// Cell center coordinates
        /// Координаты центра ячейки
        /// </summary>
        public Coordinate Coordinate { get; private set; }
        /// <summary>
        /// designated constructor
        /// </summary>
        /// <param name="coord"> Coordinate pair</param>
        /// <param name="vX"> U component</param>
        /// <param name="vY"> V component</param>
        public PrognosticCell(Coordinate coord, double vX, double vY)
        {
            this.Coordinate = coord;
            this.VelocityX = vX;
            this.VelocityY = vY;
            if (!(this.CheckValue()))
            {
                throw new System.ArgumentException("Prognostic wind velocities are incorrect");
            }
        }
        /// <summary>
        /// Проверка полей на валидность
        /// </summary>
        /// <returns></returns>
        public bool CheckValue()
        {
            bool res1 = Double.IsNaN(VelocityX);
            bool res2 = Double.IsNaN(VelocityY);
            if (res1 != res2) {
                return false;
            }
            // пустые данные
            if (res1) {
                return true;
            }
            ///Скорость ветра на высоте 10м от поверхности на Земле не может превышать 70м/c
            return Math.Abs(Math.Pow(VelocityX, 2) + Math.Pow(VelocityY, 2)) <= 4900;
        }
    }
    /// <summary>
    /// Cell obj for climate wind regular data
    /// Объект - ячейка полей ветра определенной повторяемости на 10м на регулярной сетке
    /// </summary>
    public class ClimateCell
    {
        /// <summary>
        /// once-in-5-years frequency wind, m/s 
        /// скорость ветра повторяемости один раз в 5 лет, м/с
        /// </summary>
        public double Wind5;
        /// <summary>
        /// once-in-10-years frequency wind, m/s
        /// скорость ветра повторяемости один раз в 10 лет, м/с
        /// </summary>
        public double Wind10 { get; private set; }
        /// <summary>
        /// once-in-15-years frequency wind, m/s
        /// скорость ветра повторяемости один раз в 15 лет, м/с
        /// </summary>
        public double Wind15 { get; private set; }
        /// <summary>
        /// Cell center coordinate pair 
        /// Координаты центра ячейки
        /// </summary>
        public Coordinate Coordinate { get; private set; }
        /// <summary>
        /// designated constructor
        /// </summary>
        /// <param name="coord"></param>
        /// <param name="w5"></param>
        /// <param name="w10"></param>
        /// <param name="w15"></param>
        public ClimateCell(Coordinate coord, double w5, double w10, double w15)
        {
            this.Coordinate = coord;
            this.Wind5 = w5;
            this.Wind10 = w10;
            this.Wind15 = w15;
            if (!(this.CheckValue()))
            {
                throw new System.ArgumentException("Climate wind value is not correct");
            }
        }
        /// <summary>
        /// Провекра валидности полей класса
        /// </summary>
        /// <returns></returns>
        public bool CheckValue()
        {
            bool c5 = Double.IsNaN(Wind5);
            bool c10 = Double.IsNaN(Wind10);
            bool c15 = Double.IsNaN(Wind15);
            if (c5 != c10 || c5 != c15 || c5 != c10) {
                return false;
            }
            if (c5) {
                return true;
            }
            // ветер по модулю не должен превышать 70м/с
            return Math.Abs(Wind5) < 70 && Math.Abs(Wind10) < 70 && Math.Abs(Wind15) < 70;
        }
    }
    /// <summary>
    /// Cell Size parameters
    /// Параметры ячейки (регулярной сетки)
    /// </summary>
    public class CellSize
    {
        /// <summary>
        /// ширина ячейки (расстояние между соседними по широте центрами ячеек)
        /// </summary>
        public double Width { get; set; }
        /// <summary>
        /// высота ячейки (расстояние между соседними по долготе центрами ячеек)
        /// </summary>
        public double Height { get; set; }
        /// <summary>
        /// designated constructor
        /// </summary>
        /// <param name="wdh">Cell Width</param>
        /// <param name="hgh">Cell Height</param>
        public CellSize(double wdh, double hgh)
        {
            this.Width = wdh;
            this.Height = hgh;
            if (!(this.CheckValue()))
            {
                throw new System.ArgumentException("Cell width or height values are incorrect!");
            }
        }
        /// <summary>
        /// Проверка валидности полей
        /// </summary>
        /// <returns></returns>
        public bool CheckValue()
        {
            return Width > 0 && Height > 0;
        }
    }
    /// <summary>
    /// power line object
    /// объект - ЛЭП
    /// </summary>
    public class Powerline
    {
        /// <summary>
        ///  unique id
        ///  уникальный идентификатор 
        /// </summary>
        public int Identifier { get; set; }
        /// <summary>
        /// year of construction
        /// год постройки ЛЭП
        /// </summary>
        public int Year { get; set; }
        /// <summary>
        /// average height of cable span between two poles, meters
        /// средняя высота пролета, м
        /// </summary>
        public double Height { get; set; }
        /// <summary>
        /// power kV for switches
        /// Напряжение ЛЭП, кВ
        /// </summary>
        public int Power { get; set; }
        /// <summary>
        /// Line vertices coordinate list
        /// список координат вершин ЛЭП как линейного объекта
        /// </summary>
        public IList<Coordinate> Coordinates { get; set; }
        /// <summary>
        /// assigned powerstation/pole
        /// идентификатор соответсвующего конца/начала линии (столб, трансформаторная подстанция, понижающая подстанция)
        /// </summary>
        public int PointFromID { get; set; }
        /// <summary>
        /// assigned powerstation/pole
        /// идентификатор соответсвующего конца/начала линии (столб, трансформаторная подстанция, понижающая подстанция)
        /// </summary>
        public int PointToID { get; set; }
        /// <summary>
        /// broken/not broken switch
        /// сломана (true) или нет (false) линяя
        /// </summary>
        public bool IsBroken { get; set; }
        /// <summary>
        /// power on switch
        /// получает (true) или нет (false) линия питание
        /// </summary>
        public bool IsON { get; set; }

        public Powerline()
        {
            //default constructor body
        }
        /// <summary>
        /// designated constructor
        /// </summary>
        /// <param name="coord"></param>
        /// <param name="id"></param>
        /// <param name="yer"></param>
        /// <param name="h"></param>
        /// <param name="pw"></param>
        /// <param name="isbrkn"></param>
        /// <param name="ison"></param>
        /// <param name="toID"></param>
        /// <param name="fromID"></param>
        public Powerline(IList<Coordinate> coordinates, int id, int year, double height, int power, int toID, int fromID)
        {
            this.Coordinates = coordinates;
            this.Identifier = id;
            this.Year = year;
            this.Height = height;
            this.Power = power;
            this.IsBroken = false;
            this.IsON = false;
            this.PointFromID = fromID;
            this.PointToID = toID;
            if (!(this.CheckValue()))
            { 
                throw new System.ArgumentException("Powerline object wasn't initialized correctly");
            }
        }
        /// <summary>
        /// проверка валидности полей
        /// </summary>
        /// <returns></returns>
        public bool CheckValue()
        {
            bool checker = 
                Identifier >= 0 && Math.Abs(Year - 1985) < 45 && 
                Math.Abs(Height - 15) < 15 && Power > 0 && Power < 1000 && 
                PointFromID >= 0 && PointToID >= 0;
            return checker;
        }
    }
    /// <summary>
    /// powerstation/pole point class
    /// класс для трансформаторных подстанций/столбов/понижающих(конечных) подстанций
    /// </summary>
    public class PowerStation
    {
        /// <summary>
        /// unique id  
        /// уникальный идентификатор подстанции/столба
        /// </summary>
        public int Identifier { get; set; }
        /// <summary>
        /// Coordinates
        /// </summary>
        public Coordinate Coordinate { get; set; }
        /// <summary>
        /// station name field
        /// название подстанции
        /// </summary>
        public string Name { get; set; }
        /// <summary>
        /// power, kW
        /// мощность, кВ
        /// </summary>
        public int Power { get; set; }
        /// <summary>
        /// type of point - trans/pole/endpoint
        /// тип станции - трансформаторная подстанция/столб/понижающая(конечная)подстанция
        /// </summary>
        public enum StationType 
        {
           Pole, Trans, Endstat
        };
        /// <summary>
        /// тип подстанции
        /// </summary>
        public StationType Stationtype { get; set; }
        /// <summary>
        /// is point a source?
        /// является ли подстаниция источником питания (в случае ТЭЦ или питания от внешних для рассматриваемой цепи ЛЭП)
        /// </summary>
        public bool IsSource { get; set; }
        /// <summary>
        /// power on switch
        /// поступает (true) или нет (false) на подстанцию питание
        /// </summary>
        public bool IsON { get; set; }
        /// <summary>
        /// asigned powerlines list
        /// список оканчивающихся/начинающихся на подстанции ЛЭП
        /// </summary>
        public IList<Powerline> LineList { get; set; }
        public PowerStation()
        {
            //default constructor
        }
        /// <summary>
        /// designated constructor
        /// </summary>
        /// <param name="crds"></param>
        /// <param name="id"></param>
        /// <param name="stname"></param>
        /// <param name="stpower"></param>
        /// <param name="sttype"></param>
        /// <param name="issource"></param>
        /// <param name="ison"></param>
        public PowerStation(Coordinate crds, int id, string stname, int stpower, StationType sttype, bool issource)
        {
            this.Coordinate = crds;
            this.Identifier = id;
            this.Name = stname;
            this.Power = stpower;
            this.Stationtype = sttype;
            this.IsSource = issource;
            this.IsON = false;
        }
        public bool CheckValue()
        {
            bool checker = Identifier >=0 && Power >0 && Power < 1000;
            return checker;
        }
    }

    enum FunctionType
    {
        FunctionVelocityX = 0,
        FunctionVelocityY = 1,
        FunctionClimate5 = 2,
        FunctionClimate10 = 3,
        FunctionClimate15 = 4
    }
    /// <summary>
    /// DTO for input
    /// </summary>
    public class Input
    {
        /// <summary>
        /// prognistic raster info
        /// массив прогностического ветра
        /// </summary>
        public List<List<PrognosticCell>> PrognosticCells { get; set; }
        /// <summary>
        /// prognostic raster cell info
        /// параметры ячеек регулярной сетки прогностического ветра
        /// </summary>
        public CellSize PrognosticCellSize { get; set; }
        /// <summary>
        /// affine coefficients from prognostic raster projections
        /// коэффициенты аффиного проеобразования из проекции массива прогностического ветра
        /// </summary>
        public double[] PrognosticAffineCoefficients { get; set; }

        /// <summary>
        /// climate raster array
        /// массив климатических полей скорости ветра заданной повторяемости
        /// </summary>
        public List<List<ClimateCell>> ClimateCells { get; set; }
        /// <summary>
        /// climate raster cell info
        /// параметры ячеек регулярной сетки климатических полей скорости ветра заданной повторяемости 
        /// </summary>
        public CellSize ClimateCellSize { get; set; }
        /// <summary>
        /// affine coefficients from climate raster projection
        /// коэффициенты аффинного преобразования из проекции массива климатических полей скорости ветра заданной повторяемости
        /// </summary>
        public double[] ClimateAffineCoefficients { get; set; }

        /// <summary>
        /// lines list
        /// список ЛЭП
        /// </summary>
        public List<Powerline> PowerLines { get; set; }
        /// <summary>
        /// stations/poles list
        /// список точечных объектов  - трансформаторных подстанций/столбов/понижающих(конечных) подстанций 
        /// </summary>
        public List<PowerStation> PowerStations { get; set; }
        /// <summary>
        /// maximum distance for line segment, meters
        /// максимальное расстояние между точками ЛЭП, для которых проверяется сломаются/несломаются под действием ветра
        /// </summary>
        public double DistThreshold
        {
            get
            {
                return 500;
            }
        }
    }
    /// <summary>
    /// Output
    /// </summary>
    public class Output
    {
        /// <summary>
        /// stations list without power
        /// Список подстанций на которые не поступит питание в результате предсказанных поломок ЛЭП в сети
        /// </summary>
        public List<PowerStation> DisabledStations { get; set; }
        /// <summary>
        /// broken lines list
        /// Список прогнозируемых сломанных ЛЭП в результате ветрового воздействия
        /// </summary>
        public List<Powerline> DisabledLines { get; set; }
    }
    /// <summary>
    /// main calculations class
    /// </summary>
    public class StressPowerChecker
    {
        /// <summary>
        /// Input Data
        /// </summary>
        private Input Input;

        /// <summary>
        /// Main function for power graph algorithm
        /// Основая функция для вычисления распределения питания по графу ЛЭП-подстанции
        /// </summary>
        public Output CheckPower(Input input)
        {
            this.Input = input;
            //Calculate which lines are broken
            List<WindStressPRM.Powerline> prmBrokenLines = brokenPowerLinesAfterCheck();
            //get the graph
            PreparingPowerItems();
            //start from source points
            foreach (PowerStation pwstation in input.PowerStations)
            {
                if (pwstation.IsSource)
                {
                    CheckPowerPointsForStation(pwstation);
                }
            }
            //fill output
            Output output = new Output();
            output.DisabledStations = new List<PowerStation>();
            output.DisabledLines = new List<Powerline>();
            foreach (Powerline line in input.PowerLines)
            {
                if (line.IsBroken)
                {
                    output.DisabledLines.Add(line);
                }
            }
            foreach (PowerStation powerStation in input.PowerStations)
            {
                //stations of type pole can be disabled if the line is broken 
                if (!powerStation.IsON && !(powerStation.Stationtype == PowerStation.StationType.Pole))
                {
                    output.DisabledStations.Add(powerStation);
                }
            }
            return output;
        }
        /// <summary>
        /// recursive search function for power graph
        /// </summary>
        /// <param name="sourcepoint"> powered station to search next powered station from</param>
        private void CheckPowerPointsForStation(PowerStation sourcepoint)
        {
            if (!sourcepoint.IsON)
            {
                throw new Exception("CheckPowerPointsForStation is called from disabled sourcepoint");
            }
            // if the point is not a pole  - i.e. we know 
            // it can redistribute power within connected lines
            // we turn it ON if any of the  connected lines are powered
            foreach (Powerline line in sourcepoint.LineList)
            {
                if (!line.IsBroken && !line.IsON)
                {
                    line.IsON = true;
                    foreach (PowerStation powerStation in Input.PowerStations)
                    {
                        if (powerStation.Identifier != sourcepoint.Identifier && (powerStation.Identifier == line.PointFromID || powerStation.Identifier == line.PointToID))
                        {
                            if (!(sourcepoint.Stationtype == PowerStation.StationType.Pole)) {
                                powerStation.IsON = true;
                                CheckPowerPointsForStation(powerStation);
                            }
                            else {
                                // if line is a pole we have to check if it's actually able to
                                // get electricity to other points i.e. no connected lines are broken
                                bool powerLineCheck = false;
                                foreach (Powerline powerline in powerStation.LineList)
                                {
                                    if (powerline.IsBroken) {
                                        powerLineCheck = true; 
                                    }
                                }
                                if (!powerLineCheck)
                                {
                                    powerStation.IsON = true;
                                    CheckPowerPointsForStation(powerStation);
                                }
                            }
                        }
                    }
                }
                //else we have broken line or already switched ON powerpoint
            }

        }
        /// <summary>
        /// preparing powerstations and lists of lines for them to get a graph-like structure
        /// </summary>
        private void PreparingPowerItems()
        {
            //First we make sure that all the sources are ON
            //and all non sources are OFF
            foreach (PowerStation powerStation in Input.PowerStations)
            {
                if (powerStation.IsSource == true) {
                    powerStation.IsON = true; 
                }
                else {
                    powerStation.IsON = false; 
                } 
                // for each power station we create a list of powerlines it is attached to
                List<Powerline> lines = new List<Powerline>();

                foreach (Powerline line in Input.PowerLines)
                {
                    //we also switch OFF all lines
                    line.IsON = false;
                    if (line.PointFromID == powerStation.Identifier || line.PointToID == powerStation.Identifier) {
                        lines.Add(line);
                    }
                }
                powerStation.LineList = lines;
            }
        }
        /// <summary>
        /// function to chec if powerline is broken by wind
        /// </summary>
        /// <returns>list of broken powerlines</returns>
        private List<Powerline> brokenPowerLinesAfterCheck()
        {
            List<Powerline> brokenLines = new List<Powerline>();
            // actually there are curves in powerLines
            foreach (Powerline powerCurve in Input.PowerLines)
            {
                // get coordinates list
                IList<Coordinate> points = powerCurve.Coordinates;

                List<bool> checkList = new List<bool>();
                // cycle throw all points in line
                for (int i = 1; i < points.Count; i++)
                {
                    double x1 = points[i - 1].X;
                    double y1 = points[i - 1].Y;
                    double x2 = points[i].X;
                    double y2 = points[i].Y;
                    bool result = linearLineIsBroken(points[i - 1], points[i], powerCurve.Height, powerCurve.Power);
                    checkList.Add(result);
                }
                foreach (bool chkpnt in checkList)
                {
                    if (chkpnt == true)
                    {
                        powerCurve.IsBroken = true;
                    }
                }
                if (powerCurve.IsBroken)
                {
                    brokenLines.Add(powerCurve);
                }
            }

            return brokenLines;
        }

        /// <summary>
        /// Divides the line between two vertices if needed
        /// depending on dist_threshold 
        /// and checks if line is broken by wind in any of midpoints
        /// </summary>
        /// <param name="coord1">first vertice coordinate</param>
        /// <param name="coord2">second vertice </param>
        /// <param name="heightLine">height</param>
        /// <param name="power">power for climatology switches</param>
        /// <returns>true if line is broken false otherwise</returns>
        private bool linearLineIsBroken(Coordinate coord1, Coordinate coord2, double heightLine, int power)
        {
            double distance = Math.Sqrt((coord2.X - coord1.X) * (coord2.X - coord1.X) + (coord2.Y - coord1.Y) * (coord2.Y - coord1.Y));
            double distpropD = distance / Input.DistThreshold;
            List<Coordinate> pointlist = new List<Coordinate>();
            Coordinate midpoint = new Coordinate(0, 0);
            int distpropI = Convert.ToInt32(distpropD);
            if (distpropI > 1)
            {
                double constXdiff = (coord2.X - coord1.X) / distpropI;
                double constYdiff = (coord2.Y - coord1.Y) / distpropI;
                Coordinate subCoord1 = new Coordinate(coord1.X, coord1.Y);
                Coordinate subCoord2 = new Coordinate(coord1.X + constXdiff, coord1.Y + constXdiff);
                for (int j = 0; j < distpropI; j++)
                {
                    if (j > 0)
                    {
                        // move to next center, at the next interval
                        subCoord1 = new Coordinate(subCoord2.X, subCoord2.Y);
                        subCoord2 = new Coordinate(subCoord1.X + constXdiff, subCoord1.Y + constYdiff);
                    }
                    midpoint = new Coordinate((subCoord1.X + subCoord2.X) / 2, (subCoord1.Y + subCoord2.Y) / 2);
                    pointlist.Add(midpoint);
                }
            }
            else
            {
                midpoint = new Coordinate( (coord1.X + coord2.X) / 2, (coord1.Y + coord2.Y) / 2);
                pointlist.Add(midpoint);
            }
            FunctionType climateType;
            if (power > 5 && power < 330)
            {
                climateType = FunctionType.FunctionClimate10;
            }
            else
            {
                if (power <= 5) { 
                    climateType = FunctionType.FunctionClimate5; 
                }
                else { 
                    climateType = FunctionType.FunctionClimate15;
                }
            }
            List<bool> checkbool = new List<bool>();
            foreach (Coordinate coords in pointlist)
            {
                bool res = false;
                double uwind = interpol(coords, FunctionType.FunctionVelocityX);
                double vwind = interpol(coords, FunctionType.FunctionVelocityY);
                double climwind = interpol(coords, climateType);
                double umod = Math.Sqrt(uwind * uwind + vwind * vwind); ;
                double angleline = Math.Atan2((coord2.Y - coord1.Y), (coord2.X - coord1.X));
                double anglewind = Math.Atan2(vwind, uwind) - angleline;
                double sinwind = Math.Sin(anglewind);
                double C_height = 1.0;
                if (umod < 20)
                { //wind is too low 
                    res = false;
                }
                else
                {   // calculate prognostic and threshold windstress
                    double p1 = -0.00078501;
                    double p2 = 0.13431;
                    double p3 = -2.11112;
                    double p4 = 19.548;
                    double qpr = p1 * umod * umod * umod + p2 * umod * umod + p3 * umod + p4;
                    double qcl = p1 * climwind * climwind * climwind + p2 * climwind * climwind + p3 * climwind + p4;
                    double Ppr = qpr * C_height * sinwind * sinwind;
                    double Pcl = qcl * C_height * 1.0;
                    if (Ppr >= Pcl)
                    {
                        // here line is broken
                        res = true;
                    }
                }
                checkbool.Add(res);
            }
            bool result = false;
            foreach (bool lastcheck in checkbool)
            {
                if (lastcheck == true) 
                { 
                    result = true;
                }
            }
            return result;
        }

        private double[] affineCoefficients(FunctionType functionType)
        {
            switch (functionType)
            {
                case FunctionType.FunctionVelocityX:
                case FunctionType.FunctionVelocityY:
                    {
                        return Input.PrognosticAffineCoefficients;
                    }
                case FunctionType.FunctionClimate5:
                case FunctionType.FunctionClimate10:
                case FunctionType.FunctionClimate15:
                    {
                        return Input.ClimateAffineCoefficients;
                    }
                default:
                    break;
            }
            return null;
        }

        // Taken from DotSpatial https://github.com/ViceIce/DotSpatial/blob/22c156c7646b1595d88d2523c066a9c6ab4d3a53/DotSpatial.Data/RasterBoundsExt.cs
        // RasterBoundsExt.cs. Define AffineCoefficients like this.
        private Index projectionToCell(Coordinate coordinate, FunctionType functionType)
        {
            double[] c = affineCoefficients(functionType);
            double rw, cl;
            if (c[2] == 0 && c[4] == 0)
            {
                // If we don't have skew terms, then the x and y terms are independent.
                // also both of the other tests will have divide by zero errors and fail.
                cl = (coordinate.X - c[0]) / c[1];
                rw = (coordinate.Y - c[3]) / c[5];
            }
            else if (c[2] != 0)
            {
                // this equation will have divide by zero if c[2] is zero, but works if c[4] is zero
                double div = (c[5] * c[1]) / c[2] - c[4];
                cl = (c[3] + (c[5] * coordinate.X) / c[2] - (c[5] * c[0]) / c[2] - coordinate.Y) / div;
                rw = (coordinate.X - c[0] - c[1] * cl) / c[2];
            }
            else
            {
                double div = (c[4] * c[2] / c[1] - c[5]);
                rw = (c[3] + c[4] * coordinate.X / c[1] - c[4] * c[0] / c[1] - coordinate.Y) / div;
                cl = (coordinate.X - c[2] * rw - c[0]) / c[1];
            }
            int iRow = (int)Math.Round(rw);
            int iCol = (int)Math.Round(cl);

            if (iRow < 0 || iCol < 0 || iRow >= countInList(functionType, true) || iCol >= countInList(functionType, false))
            {
                throw new Exception("projectionToCell method trying to find uncorrect index");
            }
            return new Index(iRow, iCol);
        }
        private Coordinate cellToProjection(Index index, FunctionType functionType)
        {
            switch (functionType)
            {
                case FunctionType.FunctionVelocityX:
                case FunctionType.FunctionVelocityY:
                    {
                        return Input.PrognosticCells[index.Row][index.Col].Coordinate;
                    }
                case FunctionType.FunctionClimate5:
                case FunctionType.FunctionClimate10:
                case FunctionType.FunctionClimate15:
                    {
                        return Input.ClimateCells[index.Row][index.Col].Coordinate;
                    }
                default:
                    break;
            }
            return null;
        }

        int countInList(FunctionType functionType, bool forRows)
        {
            switch (functionType)
            {
                case FunctionType.FunctionVelocityX:
                case FunctionType.FunctionVelocityY:
                    {
                        return forRows ? Input.PrognosticCells.Count : Input.PrognosticCells[0].Count;
                    }
                case FunctionType.FunctionClimate5:
                case FunctionType.FunctionClimate10:
                case FunctionType.FunctionClimate15:
                    {
                        return forRows ? Input.ClimateCells.Count : Input.ClimateCells[0].Count;
                    }
                default:
                    break;
            }
            return 0;
        }

        private double valueForFunction(FunctionType functionType, Index index)
        {
            switch (functionType)
            {
                case FunctionType.FunctionVelocityX:
                    {
                        return Input.PrognosticCells[index.Row][index.Col].VelocityX;
                    }
                case FunctionType.FunctionVelocityY:
                    {
                        return Input.PrognosticCells[index.Row][index.Col].VelocityY;
                    }
                case FunctionType.FunctionClimate5:
                    {
                        return Input.ClimateCells[index.Row][index.Col].Wind5;
                    }
                case FunctionType.FunctionClimate10:
                    {
                        return Input.ClimateCells[index.Row][index.Col].Wind10;
                    }
                case FunctionType.FunctionClimate15:
                    {
                        return Input.ClimateCells[index.Row][index.Col].Wind15;
                    }
                default:
                    break;
            }
            return 0;
        }

        private CellSize cellSizeForFunction(FunctionType functionType)
        {
            switch (functionType)
            {
                case FunctionType.FunctionVelocityX:
                case FunctionType.FunctionVelocityY:
                    {
                        return Input.PrognosticCellSize;
                    }
                case FunctionType.FunctionClimate5:
                case FunctionType.FunctionClimate10:
                case FunctionType.FunctionClimate15:
                    {
                        return Input.ClimateCellSize;
                    }
                default:
                    break;
            }
            throw new Exception("There is no cell size");
        }

        private double interpol(Coordinate coords, FunctionType functionType)
        {
            // select directions for projections
            const bool normalX = true;// true - East, false West
            const bool normalY = false;// true - North, false South
            Index rc = projectionToCell(coords, functionType);
            Coordinate center = cellToProjection(rc, functionType);
            double xDiff = coords.X - center.X;
            double yDiff = coords.Y - center.Y;
            //calculate second index
            int row2, col2;
            if ((xDiff >= 0 && normalX) || (!normalX && xDiff < 0))
            {
                row2 = rc.Row >= countInList(functionType, true) - 1 ? rc.Row - 1 : rc.Row + 1;
            }
            else
            {
                row2 = rc.Row > 0 ? rc.Row - 1 : rc.Row + 1;
            }
            if ((yDiff >= 0 && normalY) || (!normalY && yDiff < 0))
            {
                col2 = rc.Col >= countInList(functionType, false) - 1 ? rc.Col - 1 : rc.Col + 1;
            }
            else
            {
                col2 = rc.Col > 0 ? rc.Col - 1 : rc.Col + 1;
            }
            // indexes and values at bounds
            Index rcBotLeft = new Index(Math.Min(row2, rc.Row), Math.Min(col2, rc.Col));
            Index rcBotRight = new Index(Math.Max(row2, rc.Row), Math.Min(col2, rc.Col));
            Index rcTopLeft = new Index(Math.Min(row2, rc.Row), Math.Max(col2, rc.Col));
            Index rcTopRight = new Index(Math.Max(row2, rc.Row), Math.Max(col2, rc.Col));
            double valBotLeft = valueForFunction(functionType, rcBotLeft);
            double valBotRight = valueForFunction(functionType, rcBotRight);
            double valTopLeft = valueForFunction(functionType, rcTopLeft);
            double valTopRight = valueForFunction(functionType, rcTopRight);
            Coordinate origin = cellToProjection(rcBotLeft, functionType);
            
            bool testBotLeft = Double.IsNaN(valBotLeft);
            bool testBotRight = Double.IsNaN(valBotRight);
            bool testTopLeft = Double.IsNaN(valTopLeft);
            bool testTopRight = Double.IsNaN(valTopRight);
            if (testBotLeft || testBotRight || testTopLeft || testTopRight) 
            {
                // tests indicates that value at test-cell is missed. 
                if (testTopRight && testTopLeft && testBotLeft && testBotRight) 
                {
                    throw new Exception("Current value in such a bad place, you need to fill these place with raster values");
                }
                int count = 0;
                if (testBotLeft)
                {
                    valBotLeft = 0;
                    count++;
                }
                if (testBotRight)
                {
                    valBotRight = 0;
                    count++;
                }
                if (testTopLeft)
                {
                    valTopLeft = 0;
                    count++;
                }
                if (testTopRight)
                {
                    valTopRight = 0;
                    count++;
                }
                //of course there is count not 0;
                if (count == 0) {
                    throw new Exception("Interpolation Logic Error");
                }
                double average = (valTopLeft + valTopRight + valBotLeft + valBotRight)/count;
                if (testBotLeft)
                {
                    valBotLeft = average;
                }
                if (testBotRight)
                {
                    valBotRight = average;
                }
                if (testTopLeft)
                {
                    valTopLeft = average;
                }
                if (testTopRight)
                {
                    valTopRight = average;
                }