惠更斯Rothermel模型模拟-森林火灾-计算坐标（后台）

using System;
using System.Collections.Generic;
using System.Linq;
using System.Web;
using System.Web.UI;
using System.Web.UI.WebControls;

public partial class _Default : System.Web.UI.Page
{
    //变量定义
    double Wn;   //可燃物净载量
    double v;   //理想反应速率
    double Time; //时间步长
    double Toltime; //反应时间
    double density; //可燃物密度
    double h;        //可燃物烩值
    double m;        //水分对反应速率的影响系数
    double s;        //矿物资对反应速率的影响系数
    double a;        //有效加热数
    double p;        //可燃物表面积体积比
    double q;        //林火传播通率
    double Qi;       //单位可燃物的点燃热
    double x;
    double y;
    protected void Page_Load(object sender, EventArgs e)
    {

    }

    protected void Button1_Click(object sender, EventArgs e)
    {
        x = Convert.ToDouble(TextBox1.Text);
        y= Convert.ToDouble(TextBox2.Text);
        //获取数据
        Wn = Convert.ToDouble(TextBox6.Text);
        density = Convert.ToDouble(TextBox7.Text);
        v = Convert.ToDouble(TextBox8.Text);
        Time = Convert.ToDouble(TextBox12.Text) + 10;

        Toltime = Convert.ToDouble(TextBox11.Text) - 10;
        Toltime = Toltime - 30;
        h = 15000;
        m = 0.6;
        s = 0.8;
        a = 0.6;
        p = 2;
        q = 0.8;
        Qi = 2400;
        class1 c12 = new class1();
        int num = (int)(Toltime / Time);
        for (int i = 0; i < 8; i++)
        {
            HuygenSpread(x, y, num, i, 0);
        }
        for (int i = 0; i < list1.Count; i++)
        {
            class1 c13 = new class1();
            c13 =list1[i];
            Response.Write("["+c13.X+","+ c13.Y+"],");
        }
    }
    经纬度转墨卡托
    public class1 lonLat2Mercator(double x1,double y1)
    {
        class1 c1 = new class1();
        c1.X = x1 * 20037508.34 / 180;
        double Y2 = Math.Log(Math.Tan((90 + y1) * Math.PI / 360)) / (Math.PI / 180);
        c1.Y = Y2 * 20037508.34 / 180;
        return c1;
    }
    //墨卡托转经纬度
    public class1 Mercator2lonLat(double x1, double y1)
    {
        class1 c1 = new class1();
        c1.X = x1 / 20037508.34 * 180;
        double Y2 = y1 / 20037508.34 * 180;
        c1.Y = 180 / Math.PI * (2 * Math.Atan(Math.Exp(Y2 * Math.PI / 180)) - Math.PI / 2);
        return c1;
    }
    List<class1> list1 = new List<class1>();
    List<class1> list2 = new List<class1>();
    public void HuygenSpread(double x1, double y1, int m, int n, int i)
    {
        class1 classxy = new class1();
        double s;
        double v;
        double x;
        double y;
        int j = m;
        int k = n;
        object missing = Type.Missing;
        if (i < j)
        {
            v = SpreadSpeed(k, Convert.ToDouble(TextBox4.Text));
            s = v * Time;
            //  x = (double)pFeature.get_Value(pFeature.Fields.FindField("X_"));
            // y = (double)pFeature.get_Value(pFeature.Fields.FindField("Y_"));
            class1 c2 = new class1();
            c2 = lonLat2Mercator(x1, y1);
            class1 c1 = GetNextFeature(c2.X, c2.Y, 1.4 * s, k, m);
            c1 = Mercator2lonLat(c1.X, c1.Y);
            list1.Add(c1);
        }
        else
        {
            class1 c3 = new class1();
            c3.X = x1;
            c3.Y = y1;
            list1.Add(c3);

        }
         return;
      

    
    }
    //获取坐标
    public class1 GetNextFeature(double x, double y, double s, int n,int sj)
    {
        double s1;
        int m;
        class1 c1 = new class1();
        class1 c2 = new class1();
        double fengxiang = Convert.ToDouble(TextBox4.Text);
        int number = 80;
        switch (n)
        {
            case 0:
                {
                    s1 = s * Math.Sin((Math.PI) / 4);
                    m = (int)(s1 / 80);
                    // c1.X = (x - m * 80 - 0.1 - number * sj);
                    //c1.Y = (y + m * 80 + 0.1 + number * sj);
                    if (fengxiang > 0 && fengxiang < 90)
                    {
                        c1.X = (x - m * 80 - 0.1 - number * sj);
                        c1.Y = (y + m * 80 + 0.1 + number * sj);
                    }
                    else if (fengxiang > 90 && fengxiang < 180)
                    {
                        c1.X = (x - m * 80 - 0.1 - number * sj-80*2);
                        c1.Y = (y + m * 80 + 0.1 + number * sj+80*2);
                    }
                    else if (fengxiang > 180 && fengxiang < 270)
                    {
                        c1.X = (x - m * 80 - 0.1 - number * sj - 80 * 3);
                        c1.Y = (y + m * 80 + 0.1 + number * sj + 80 * 3);
                    }
                    else if (fengxiang > 270 && fengxiang < 360)
                    {
                        c1.X = (x - m * 80 - 0.1 - number * sj-80*1);
                        c1.Y = (y + m * 80 + 0.1 + number * sj+80*1);
                    }
                    else if (fengxiang == 0|| fengxiang == 360)
                    {
                        c1.X = (x - m * 80 - 0.1 - number * sj);
                        c1.Y = (y + m * 80 + 0.1 + number * sj);
                    }
                    else if (fengxiang == 90)
                    {
                        c1.X = (x - m * 80 - 0.1 - number * sj);
                        c1.Y = (y + m * 80 + 0.1 + number * sj);
                    }
                    else if (fengxiang == 180)
                    {
                        c1.X = (x - m * 80 - 0.1 - number * sj-80*1);
                        c1.Y = (y + m * 80 + 0.1 + number * sj+80*1);
                    }
                    else if (fengxiang == 270)
                    {
                        c1.X = (x - m * 80 - 0.1 - number * sj - 80*2);
                        c1.Y = (y + m * 80 + 0.1 + number * sj + 80*2);
                    }
                    break;
                }
            case 1:
                {
                    s1 = s;
                    m = (int)(s1 / 80);
                    // c1.X = (x - 0.1 - number * sj);
                    // c1.Y = (y + m * 80 + number * sj + 0.1);
                    if (fengxiang > 0 && fengxiang < 90)
                    {
                         c1.X = (x - 0.1-80);
                         c1.Y = (y + m * 80 + number * sj + 0.1);
                    }
                    else if (fengxiang > 90 && fengxiang < 180)
                    {
                        c1.X = (x - 0.1 - 80);
                        c1.Y = (y + m * 80 + number * sj + 0.1+80*1);
                    }
                    else if (fengxiang > 180 && fengxiang < 270)
                    {
                        c1.X = (x - 0.1 - 80);
                        c1.Y = (y + m * 80 + number * sj + 0.1+ 80 * 2);
                    }
                    else if (fengxiang > 270 && fengxiang < 360)
                    {
                        c1.X = (x - 0.1 - 80);
                        c1.Y = (y + m * 80 + number * sj + 0.1 +80 * 2);
                    }
                    else if (fengxiang == 0 || fengxiang == 360)
                    {
                        c1.X = (x - 0.1 - 80);
                        c1.Y = (y + m * 80 + number * sj + 0.1+80);
                    }
                    else if (fengxiang == 90)
                    {
                        c1.X = (x - 0.1 - 80);
                        c1.Y = (y + m * 80 + number * sj + 0.1);
                    }
                    else if (fengxiang == 180)
                    {
                        c1.X = (x - 0.1 - 80);
                        c1.Y = (y + m * 80 + number * sj + 0.1);
                    }
                    else if (fengxiang == 270)
                    {
                        c1.X = (x - 0.1 - 80);
                        c1.Y = (y + m * 80 + number * sj + 0.1+80*3);
                    }
                    break;
                }
            case 2:
                {
                    s1 = s * Math.Sin((Math.PI) / 4);
                    m = (int)(s1 / 80);
                    // c1.X = (x + m * 80 + number * sj + 0.1);
                    // c1.Y = (y + m * 80 + number * sj + 0.1);
                    if (fengxiang > 0 && fengxiang < 90)
                    {
                        c1.X = (x + m * 80 + number * sj + 0.1);
                        c1.Y = (y + m * 80 + number * sj + 0.1);
                    }
                    else if (fengxiang > 90 && fengxiang < 180)
                    {
                        c1.X = (x + m * 80 + number * sj + 0.1-80*2);
                        c1.Y = (y + m * 80 + number * sj + 0.1-80*2);
                    }
                    else if (fengxiang > 180 && fengxiang < 270)
                    {
                        c1.X = (x + m * 80 + number * sj + 0.1 - 80 * 1);
                        c1.Y = (y + m * 80 + number * sj + 0.1 - 80 * 1);
                    }
                    else if (fengxiang > 270 && fengxiang < 360)
                    {
                        c1.X = (x + m * 80 + number * sj + 0.1+80*5);
                        c1.Y = (y + m * 80 + number * sj + 0.1+80*5);
                    }
                    else if (fengxiang == 0 || fengxiang == 360)
                    {
                        c1.X = (x + m * 80 + number * sj + 0.1 + 80*3);
                        c1.Y = (y + m * 80 + number * sj + 0.1 + 80*3);
                    }
                    else if (fengxiang == 90)
                    {
                        c1.X = (x + m * 80 + number * sj + 0.1 + 80);
                        c1.Y = (y + m * 80 + number * sj + 0.1 + 80);
                    }
                    else if (fengxiang == 180)
                    {
                        c1.X = (x + m * 80 + number * sj + 0.1);
                        c1.Y = (y + m * 80 + number * sj + 0.1);
                    }
                    else if (fengxiang == 270)
                    {
                        c1.X = (x + m * 80 + number * sj + 0.1 + 80*2);
                        c1.Y = (y + m * 80 + number * sj + 0.1 + 80*2);
                    }
                    break;
                }
            case 3:
                {
                    s1 = s;
                    m = (int)(s1 / 80);
                   // c1.X = (x + m * 80 + 0.1 + number * sj);
                   // c1.Y = (y + 0.1 + number * sj);
                    if (fengxiang > 0 && fengxiang < 90)
                    {
                        c1.X = (x + m * 80 + 0.1 + number * sj + 80 * 2);
                        c1.Y = (y + 0.1 +80);
                    }
                    else if (fengxiang > 90 && fengxiang < 180)
                    {
                        c1.X = (x + m * 80 + 0.1 + number * sj);
                        c1.Y = (y + 0.1 + 80);
                    }
                    else if (fengxiang > 180 && fengxiang < 270)
                    {
                        c1.X = (x + m * 80 + 0.1 + number * sj);
                        c1.Y = (y + 0.1 + 80);
                    }
                    else if (fengxiang > 270 && fengxiang < 360)
                    {
                        c1.X = (x + m * 80 + 0.1 + number * sj+80 * 2);
                        c1.Y = (y + 0.1 + 80);
                    }
                    else if (fengxiang == 0 || fengxiang == 360)
                    {
                        c1.X = (x + m * 80 + 0.1 + number * sj + 80 * 3);
                        c1.Y = (y + 0.1 + 80);
                    }
                    else if (fengxiang == 90)
                    {
                        c1.X = (x + m * 80 + 0.1 + number * sj + 80 * 1);
                        c1.Y = (y + 0.1 + 80);
                    }
                    else if (fengxiang == 180)
                    {
                        c1.X = (x + m * 80 + 0.1 + number * sj);
                        c1.Y = (y + 0.1 + 80);
                    }
                    else if (fengxiang == 270)
                    {
                        c1.X = (x + m * 80 + 0.1 + number * sj + 80 *1);
                        c1.Y = (y + 0.1 + 80);
                    }
                    break;
                }
            case 4:
                {
                    s1 = s * Math.Sin((Math.PI) / 4);
                    m = (int)((s1 / 80));
                    // c1.X = (x + m * 80 + number * sj + 0.1);
                    // c1.Y = (y - m * 80 - number * sj - 0.1);
                    if (fengxiang > 0 && fengxiang < 90)
                    {
                        c1.X = (x + m * 80 + number * sj + 80 * 3);
                        c1.Y = (y - m * 80 - number * sj - 80 * 3);
                    }
                    else if (fengxiang > 90 && fengxiang < 180)
                    {
                         c1.X = (x + m * 80 + number * sj + 0.1);
                         c1.Y = (y - m * 80 - number * sj - 0.1);
                    }
                    else if (fengxiang > 180 && fengxiang < 270)
                    {
                         c1.X = (x + m * 80 + number * sj + 0.1);
                         c1.Y = (y - m * 80 - number * sj - 0.1);
                    }
                    else if (fengxiang > 270 && fengxiang < 360)
                    {
                         c1.X = (x + m * 80 + number * sj + 0.1);
                         c1.Y = (y - m * 80 - number * sj - 0.1);
                    }
                    else if (fengxiang == 0 || fengxiang == 360)
                    {
                        c1.X = (x + m * 80 + number * sj + 80 * 2);
                        c1.Y = (y - m * 80 - number * sj - 80 * 2);
                    }
                    else if (fengxiang == 90)
                    {
                        c1.X = (x + m * 80 + number * sj + 80 * 2);
                        c1.Y = (y - m * 80 - number * sj - 80 * 2);
                    }
                    else if (fengxiang == 180)
                    {
                        c1.X = (x + m * 80 + number * sj + 0.1);
                        c1.Y = (y - m * 80 - number * sj - 0.1);
                    }
                    else if (fengxiang == 270)
                    {
                        c1.X = (x + m * 80 + number * sj + 0.1);
                        c1.Y = (y - m * 80 - number * sj - 0.1);
                    }
                    break;
                }
            case 5:
                {
                    s1 = s;
                    m = (int)(s1 / 80);
                    //c1.X = (x + 0.1 + number * sj);
                    // c1.Y = (y - m * 80 - number * sj - 0.1);   
                    if (fengxiang > 0 && fengxiang < 90)
                    {
                        c1.X = (x + 0.1 + 80);
                        c1.Y = (y - m * 80 - number * sj - 0.1- 80 * 2);
                    }
                    else if (fengxiang > 90 && fengxiang < 180)
                    {
                        c1.X = (x + 0.1 + 80);
                        c1.Y = (y - m * 80 - number * sj - 0.1 - 80 * 1);
                    }
                    else if (fengxiang > 180 && fengxiang < 270)
                    {
                       c1.X = (x + 0.1 + 80);
                       c1.Y = (y - m * 80 - number * sj - 0.1);    ;
                    }
                    else if (fengxiang > 270 && fengxiang < 360)
                    {
                        c1.X = (x + 0.1 + 80);
                        c1.Y = (y - m * 80 - number * sj - 0.1);
                    }
                    else if (fengxiang == 0 || fengxiang == 360)
                    {
                        c1.X = (x + 0.1 + 80);
                        c1.Y = (y - m * 80 - number * sj - 0.1 - 80 * 1);
                    }
                    else if (fengxiang == 90)
                    {
                        c1.X = (x + 0.1 + 80);
                        c1.Y = (y - m * 80 - number * sj - 0.1 - 80 * 3);
                    }
                    else if (fengxiang == 180)
                    {
                        c1.X = (x + 0.1 + 80);
                        c1.Y = (y - m * 80 - number * sj - 0.1 - 80 * 1);
                    }
                    else if (fengxiang == 270)
                    {
                        c1.X = (x + 0.1 + 80);
                        c1.Y = (y - m * 80 - number * sj - 0.1);
                    }
                    break;
                }
            case 6:
                {
                    s1 = s * Math.Sin((Math.PI) / 4);
                    m = (int)((s1 / 80));
                    // c1.X = (x - m * 80 - 0.1 - number * sj);
                    //c1.Y = (y - m * 80 - 0.1 - number * sj);
                    if (fengxiang > 0 && fengxiang < 90)
                    {
                        c1.X = (x - m * 80 - 0.1 - number * sj);
                        c1.Y = (y - m * 80 - 0.1 - number * sj);
                    }
                    else if (fengxiang > 90 && fengxiang < 180)
                    {
                        c1.X = (x - m * 80 - 0.1 - number * sj - 80 * 3);
                        c1.Y = (y - m * 80 - 0.1 - number * sj - 80 * 3);
                    }
                    else if (fengxiang > 180 && fengxiang < 270)
                    {
                        c1.X = (x - m * 80 - 0.1 - number * sj- 80 * 1);
                        c1.Y = (y - m * 80 - 0.1 - number * sj-80 * 1);
                    }
                    else if (fengxiang > 270 && fengxiang < 360)
                    {
                        c1.X = (x - m * 80 - 0.1 - number * sj);
                        c1.Y = (y - m * 80 - 0.1 - number * sj);
                    }
                    else if (fengxiang == 0 || fengxiang == 360)
                    {
                        c1.X = (x - m * 80 - 0.1 - number * sj);
                        c1.Y = (y - m * 80 - 0.1 - number * sj);
                    }
                    else if (fengxiang == 90)
                    {
                        c1.X = (x - m * 80 - 0.1 - number * sj - 80 * 2);
                        c1.Y = (y - m * 80 - 0.1 - number * sj - 80 * 2);
                    }
                    else if (fengxiang == 180)
                    {
                        c1.X = (x - m * 80 - 0.1 - number * sj - 80 * 1);
                        c1.Y = (y - m * 80 - 0.1 - number * sj - 80 * 1);
                    }
                    else if (fengxiang == 270)
                    {
                        c1.X = (x - m * 80 - 0.1 - number * sj);
                        c1.Y = (y - m * 80 - 0.1 - number * sj);
                    }
                    break;
                }
            case 7:
                {
                    s1 = s;
                    m = (int)(s1 / 80);
                    //c1.X = (x - m * 80 - number * sj);
                    // c1.Y = (y - 0.1 - number * sj);
                    if (fengxiang > 0 && fengxiang < 90)
                    {
                        c1.X = (x - m * 80 - number * sj);
                        c1.Y = (y - 0.1-80);
                    }
                    else if (fengxiang > 90 && fengxiang < 180)
                    {
                        c1.X = (x - m * 80 - number * sj- 80 * 3);
                        c1.Y = (y - 0.1 - 80);
                    }
                    else if (fengxiang > 180 && fengxiang < 270)
                    {
                        c1.X = (x - m * 80 - number * sj - 80 * 2);
                        c1.Y = (y - 0.1 - 80);
                    }
                    else if (fengxiang > 270 && fengxiang < 360)
                    {
                        c1.X = (x - m * 80 - number * sj);
                        c1.Y = (y - 0.1 - 80);
                    }
                    else if (fengxiang == 0 || fengxiang == 360)
                    {
                        c1.X = (x - m * 80 - number * sj);
                        c1.Y = (y - 0.1 - 80);
                    }
                    else if (fengxiang == 90)
                    {
                        c1.X = (x - m * 80 - number * sj - 80 * 1);
                        c1.Y = (y - 0.1 - 80);
                    }
                    else if (fengxiang == 180)
                    {
                        c1.X = (x - m * 80 - number * sj - 80 * 3);
                        c1.Y = (y - 0.1 - 80);
                    }
                    else if (fengxiang == 270)
                    {
                        c1.X = (x - m * 80 - number * sj- 80 * 1);
                        c1.Y = (y - 0.1 - 80);
                    }
                    break;
                }
        }

        return c1;
    }
     double K = 1;//惠更斯模型修正系数
   /// <summary>
   /// 计算风速
   /// </summary>
   /// <param name="n"></param>
   /// <param name="Winddir"></param>
   /// <returns></returns>
    public double SpreadSpeed(int n, double Winddir)
    {
        double angle = 0;

        double Ir;
        double R;
        double w;
        double s;

        if (n == 0)
        {
            angle = Math.Abs(Winddir - 225);
            return K * (16 + 16 * (Math.Cos(angle * Math.PI / 180)));
        }
        if (n == 1)
        {
            angle = Math.Abs(Winddir - 270);
            return K * (15 + 15 * (Math.Cos(angle * Math.PI / 180)));
        }
        if (n == 2)
        {
            angle = Math.Abs(Winddir - 315) + 180;
        }
        if (n == 3)
        {
            angle = Math.Abs(Winddir);
        }
        if (n == 4)
        {
            angle = Math.Abs(Winddir - 45);
        }
        if (n == 5)
        {
            angle = Math.Abs(Winddir - 90);
        }
        if (n == 6)
        {
            angle = Math.Abs(Winddir - 135);
            return K * (12 + 12 * (Math.Cos(angle * Math.PI / 180)));
        }
        if (n == 7)
        {
            angle = Math.Abs(Winddir - 180);
            return K * (13 + 13 * (Math.Cos(angle * Math.PI / 180)));
        }
        w = CalWindCorr(Convert.ToDouble(TextBox13.Text));
        s = CalSlopeCorr(Convert.ToDouble(TextBox9.Text));
        Ir = CalIr();
        R = (Ir * q * (1 + w + s)) / (density * a * Qi);
        return K * (R + R * Math.Cos(angle * Math.PI / 180));
    }
    //计算火焰反应强度
    private double CalIr()
    {
        double Ir;
        Ir = v * Wn * h * m * s;
        return Ir;
    }
    //计算风速修正系数
    private double CalWindCorr(double pFeature)
    {
        double windspeed;
        double Qw;
        double B;
        double C;
        double E;
        double b;  //燃料最佳紧密度
        double c;
        c = 1.2;
        windspeed = (double)pFeature;
        E = 0.715 * (Math.Pow(Math.E, (-3.95 * 10000 * p)));
        C = 7.47 * (Math.Pow(Math.E, (-0.133 * Math.Pow(p, 0.55))));
        B = 0.02526 * Math.Pow(p, 0.54);
        b = 3.348 * Math.Pow(p, -0.8189);
        Qw = C * Math.Pow(3.281 * windspeed, B) * Math.Pow(c / b, -E);
        return Qw;
    }
    //计算坡度修正系数
    public double CalSlopeCorr(double pFeature)
    {
        double slope;
        double Qs;
        slope = (double)pFeature;
        Qs = 5.275 * Math.Pow(p, -0.3) * (Math.Tan(slope * Math.PI / 180)) * (Math.Tan(slope * Math.PI / 180));
        return Qs;
    }
}
public  class class1
{
    double x;
    double y;

    public double X
    {
        get
        {
            return x;
        }

        set
        {
            x = value;
        }
    }

    public double Y
    {
        get
        {
            return y;
        }

        set
        {
            y = value;
        }
    }
}