/*  由经纬度坐标计算面积和长度      
      完成时间：20050727 wanggang
  */

  //*************************全局变量*************************************
  var a = 6378137;                //地球长半轴
  var f = 1/298.257223563;        //地球偏心率
  var e =  Math.sqrt(2*f-f*f);
  var PI = 3.141592653589793238;  //圆周率
  var dblDF = PI/180;             //由角度到弧度转换常数
  var	EPS = 0.000000000005;
  //**********************************************************************

  //************************Lambert投影***********************************
  function LambertPrj(center_long, center_lat,false_east,false_north)  //定义Lambert投影对象
  {
    this.Lon0=center_long;                                             //中央经线
    this.Lat0=center_lat;                                              //中央纬线
    this.false_east=false_east;                                        //横向偏移
    this.false_north=false_north;                                      //纵向偏移
  }  
  
  function GetQ(lat)
  {
		var p=lat;
		var x;
		a1=1-e*e;
		a2=1-e*e*Math.sin(p)*Math.sin(p);
		a3=1/(2*e);
		a4=1-e*Math.sin(p);
		a5=1+e*Math.sin(p);
		   
		x=a1*(Math.sin(p)/a2-a3*Math.log(a4/a5));
    return (x);
  }

  LambertPrj.prototype.GetXY=function(lon,lat)                          //由经纬度计算平面坐标
  {
     var arrayXY=new Array();
     
     var qp=GetQ(PI/2);
     var q0=GetQ(this.Lat0);
     var q=GetQ(lat);
     
     var B0=Math.asin(q0/qp);
     var Bt=Math.asin(q/qp);
     var Rq=a*Math.sqrt(qp/2);
     
     var D=a*(Math.cos(this.Lat0)/Math.sqrt(1-e*e*Math.sin(this.Lat0)*Math.sin(this.Lat0)))/(Rq*Math.cos(B0));
     var B=Rq*Math.sqrt(2/(1+Math.sin(B0)*Math.sin(Bt)+Math.cos(B0)*Math.cos(Bt)*Math.cos(lon-this.Lon0)));     
     //alert("qp="+qp+"  q0="+q0+"  q="+q+"  Rq="+Rq+"  B0="+B0+"  Bt="+Bt+"  D="+D+"  B="+B);     
     arrayXY[0]=this.false_east+(B*D)*Math.cos(Bt)*Math.sin(lon-this.Lon0);
     arrayXY[1]=this.false_north+(B/D)*(Math.cos(B0)*Math.sin(Bt)-Math.sin(B0)*Math.cos(Bt)*Math.cos(lon-this.Lon0)); 
		return arrayXY;
  }
  
  LambertPrj.prototype.GetLonAndLat=function(x,y)                    //由平面坐标计算经纬度
  {
     var arrayLonLat=new Array();
     var qp=GetQ(PI/2);
     var q0=GetQ(this.Lat0);

     var B0=Math.asin(q0/qp);
     var Rq=a*Math.sqrt(qp/2);
     var D=a*(Math.cos(this.Lat0)/Math.sqrt(1-e*e*Math.sin(this.Lat0)*Math.sin(this.Lat0)))/(Rq*Math.cos(B0));

     var P=Math.sqrt(((x-this.false_east)/D)*((x-this.false_east)/D)+(D*(y-this.false_north))*(D*(y-this.false_north)));
		 var C=2*Math.asin(P/(2*Rq));
		 var B1=Math.asin(Math.cos(C)*Math.sin(B0)+(D*(y-this.false_north)*Math.sin(C)*Math.cos(B0))/P); 	 
		 arrayLonLat[0]=this.Lon0+Math.atan((x-this.false_east)*Math.sin(C)/(D*P*Math.cos(B0)*Math.cos(C)-D*D*(y-this.false_north)*Math.sin(B0)*Math.sin(C)));      
     arrayLonLat[1]=B1+(Math.pow(e,2)/3+31*Math.pow(e,4)/180+517*Math.pow(e,6)/5040)*Math.sin(2*B1)+(23*Math.pow(e,4)/360+251*Math.pow(e,6)/3780)*Math.sin(4*B1)+(761*Math.pow(e,6)/45360)*Math.sin(6*B1);
     return arrayLonLat;
  } 

  //定义该点所在的地图投影：中央经线，中央纬线，横向偏移，纵向偏移  
  var NowLambert=new LambertPrj(110 * dblDF,21 * dblDF,0,0);
	//***********************************************************************************
  
	//************************************计算面积***************************************
  //由经纬度坐标数组计算面积
  function GetAreaFromLatAndLon(aryGBPoints)
  {
    var aryPoints=ConvertGBPtsToJavaPts(aryGBPoints)
    if (aryPoints.length < 3) return -1;
	
	var NewPts = ConvertPtsCoods(aryPoints);
    var lngSumPts = NewPts.length;
		var aryFromPt = new Array();
		var aryToPt = new Array();
		var dblResultArea=0;

		for (var i=0;i<lngSumPts;i++)
		{		
			aryFromPt=NewPts[i];

			if (i==lngSumPts - 1)
				aryToPt=NewPts[0];
			else
				aryToPt=NewPts[i+1];

			dblResultArea += ((aryFromPt[1] + aryToPt[1]) * (aryToPt[0] - aryFromPt[0]))/2 ;
		}

		return Math.abs(dblResultArea);
  }
  
	 //将某经纬度坐标数组转换成平面坐标
	function ConvertPtsCoods(aryLatLonPts)
	{
		var lngSumPts = aryLatLonPts.length;
		var aryPts = new Array();
  
	  for (var i=0; i<lngSumPts; i++)
		{
			aryPts[i] = NowLambert.GetXY(aryLatLonPts[i][0]* dblDF ,aryLatLonPts[i][1] * dblDF);
		}
   
		return aryPts;
	}
  //***************************************************************************************

  //**************************************计算长度*****************************************
	//由经纬度坐标数组计算长度
  function GetLengthFromLatAndLon(aryGBPoints)
  {
   	var aryPoints=ConvertGBPtsToJavaPts(aryGBPoints)
    var lngSumPts = aryPoints.length;
		var dblResultLength=0;

    if (lngSumPts < 2) return -1;

		for (var i=0;i<lngSumPts-1;i++)
		{				
			dblResultLength += GetLengthFromTwoPt(aryPoints[i][1]*dblDF,aryPoints[i][0]*dblDF,aryPoints[i+1][1]*dblDF,aryPoints[i+1][0]*dblDF);
		}

		return dblResultLength;
  }

	//计算两个经纬度坐标点之间的长度
	 function GetLengthFromTwoPt(lat1,lon1,lat2,lon2)
	 {
     var distance = 0.0;
     var faz, baz;
     var r = 1.0 - f;
     var tu1, tu2, cu1, su1, cu2, x, sx, cx, sy, cy, y, sa, c2a, cz, e, c, d;
     var cosy1, cosy2;

		 if ( (lon1 == lon2) && (lat1 == lat2)) {
			 return distance;
		 }

     cosy1 = Math.cos(lat1);
     cosy2 = Math.cos(lat2);

     if (cosy1 == 0.0) {
       cosy1 = 0.0000000001;
     }

     if (cosy2 == 0.0) {
       cosy2 = 0.0000000001;
     }

     tu1 = r * Math.sin(lat1) / cosy1;
     tu2 = r * Math.sin(lat2) / cosy2;
     cu1 = 1.0 / Math.sqrt(tu1 * tu1 + 1.0);
     su1 = cu1 * tu1;
     cu2 = 1.0 / Math.sqrt(tu2 * tu2 + 1.0);
     x = lon2 - lon1;

     distance = cu1 * cu2;
     baz = distance * tu2;
     faz = baz * tu1;
 
		 do {
       sx = Math.sin(x);
       cx = Math.cos(x);
       tu1 = cu2 * sx;
       tu2 = baz - su1 * cu2 * cx;
       sy = Math.sqrt(tu1 * tu1 + tu2 * tu2);
       cy = distance * cx + faz;
       y = Math.atan2(sy, cy);
       sa = distance * sx / sy;
       c2a = -sa * sa + 1.0;
       cz = faz + faz;
       if (c2a > 0.0) {
         cz = -cz / c2a + cy;
       }
       e = cz * cz * 2. - 1.0;
       c = ( ( -3.0 * c2a + 4.0) * f + 4.0) * c2a * f / 16.0;
       d = x;
       x = ( (e * cy * c + cz) * sy * c + y) * sa;
       x = (1.0 - c) * x * f + lon2 - lon1;
     }
     while (Math.abs(d - x) > EPS);

     x = Math.sqrt( (1.0 / r / r - 1.0) * c2a + 1.0) + 1.0;
     x = (x - 2.0) / x;
     c = 1.0 - x;
     c = (x * x / 4.0 + 1.0) / c;
     d = (0.375 * x * x - 1.0) * x;
     x = e * cy;
     distance = 1.0 - e - e;
     
		 distance = ( ( ( (sy * sy * 4.0 - 3.0) * distance * cz * d / 6.0 - x) * d / 4.0 + cz) * sy * d + y) * c * a * r;

     return distance;
  }
  //********************************************************************************************
 
  //**************************将GeoBeans数组转换成javascript数组********************************
	function ConvertGBPtsToJavaPts(pPoints)
	{	   
	   var pNewPts=new Array();
	   var nCount=pPoints.size();
		
	   for(var i=0; i<nCount; i++)
	   {
		  	var pPoint=pPoints.dotAt(i);

		    pNewPts[i]= new Array();
			pNewPts[i][0]=pPoint.m_x;
			pNewPts[i][1]=pPoint.m_y;	
	   }

       return pNewPts;
	}  
  //********************************************************************************************