//################################################
function ellRadien(phi) {
//----------------------------------
// Datum WGS84
A  = 6378137 ;
B  = 6356752.3142 ;
e2 = 0.00669437999014 ; 

len = phi.length ;
phim = 0;
//------------------------
for (i=0; i<len; i++) {   
   
   phim += phi[i];
}  
phim = phim/len ; 
//-----------------------------------
// Krümmungsradien
   
nenner = Math.sqrt( 1.0 - e2 *( Math.pow(Math.sin(phim),2) ) ) ;
M  = A*(1.0-e2)/(Math.pow(nenner,3)) ;
N  = A/nenner ;
Rg = Math.sqrt(M*N) ;

R = Array(M,N,Rg);
return R ;
}
//################################################
function EllDistance(phi1,lam1,phi2,lam2) {

phi1 = phi1 * Math.PI/180 ;
lam1 = lam1 * Math.PI/180 ;
phi2 = phi2 * Math.PI/180 ;
lam2 = lam2 * Math.PI/180 ;
//----------------------------------
// Mittlere Breiten
//----------------------------------
phim = (phi1 + phi2) /2 ;
dphi = phi2 - phi1 ;
dlam = lam2 - lam1 ;

//----------------------------------
// Datum WGS84
A  = 6378137 ;
B  = 6356752.3142 ;
e2 = 0.00669437999014 ;   
//-----------------------------------
// Krümmungsradien
   
nenner = Math.sqrt( 1.0 - e2 *( Math.pow(Math.sin(phim),2) ) ) ;
M  = A*(1.0-e2)/(Math.pow(nenner,3)) ;
N  = A/nenner ;
Rg = Math.sqrt(M*N) ;


vq = N/N;
eq = vq - 1.0 ;
  
sphim = Math.sin(phim) ;
cphim = Math.cos(phim) ;
tphim = Math.tan(phim) ;
 //------------------------------------
	
ssin1 = ( (1-9.0*eq*Math.pow(tphim,2))/(24.0*1) )* Math.pow(dphi,2) ;
ssin2 = 1.0 - Math.pow(dlam*sphim,2)/ 24.0 ;
ssin  = N* dlam*cphim*(ssin1+ssin2) ;

//-----------------------------------
scos1 = eq * (Math.pow(tphim,2) - 1.0)*Math.pow(phim,2)/(8.0) ;
scos2 = (2.0*1+3.0 * Math.pow(tphim,2)) * Math.pow(dlam*cphim,2)/ 24.0 ;
scos  = M * dphi * (1.0-scos2-scos1) ;
   //-------------------------------------
da = (3.0+8.0*eq)* Math.pow(dphi,2)/(24.0) ;
da = da + ( Math.pow(dlam*cphim,2) / 12.0) + 1 ;
da = da* dlam*sphim ;
   //---------------------------------------
   //ix = find(scos);
   //az = ones(size(scos))*pi/2 ;
   
az = Math.atan (ssin/scos) ;
   //if scos ==0,
   //    az = pi/2 ;
   //else
   //    az = atan (ssin./scos) ;
   //end
   
if (ssin < 0.0){
	az = az + 2*Math.pi ;
}
if (scos < 0.0) {
	az = az + Math.pi ;
}
if (ssin < 0.0 && scos < 0.0) {
		az = az - 2*Math.pi ;
}
//----------------------------------------------
// Strecke s12
s12 = Math.sqrt(Math.pow(ssin,2) + Math.pow(scos,2)) ;
//--------------------------------------------  
	a12 = az - da/2.0 ;
	a21 = az + da/2.0 + Math.pi ;
 
    if (a21 > 2*Math.pi) {
        a21 = a21-2*Math.pi;
    }
	//-------------------------
    // Umwandlung in Altgrad
    //a21 = a21 *180/pi ;
    //a12 = a12 *180/pi ;
    
    return s12;
 }

 
  //####################################################
 
function flaeche(lon, lat, lonm, latm) {

   len = lon.length ;

   
  
   dKM = 111.182 ;
   var F = 0 ;
   var factor = 1000 ;
   //------------------------
   for (i=0; i<len-1; i++) {
   	
     /*P1 = Ell2Gau (lon[i]-lonm,lat[i]) ;
     P2 = Ell2Gau (lon[i+1]-lonm,lat[i+1]) ; 	
     */
     x1 = (lon[i]-lonm)   * Math.cos(lat[i]  *Math.PI/180 ) * dKM ;
     x2 = (lon[i+1]-lonm) * Math.cos(lat[i+1]*Math.PI/180 ) * dKM ;
     y1 = (lat[i]-latm) * dKM ;
     y2 = (lat[i+1]-latm) * dKM ;
     /*x1 = P1[0] /factor  ;
     y1 = P1[1] /factor ;
     x2 = P2[0] /factor ;
     y2 = P2[1] /factor ;
      */
      F += x1 * y2 - x2 * y1;
      
   }  
   F = Math.abs(F)/2 ;
   
   return F;
}
//##################################################################

//     Meridianbogenlänge

// Reihenentwicklung nach F.R. Helmert ((Großmann S.18)
function MerBogLaenge (phi) {

phi = phi * Math.PI/180 ;
//----------------------------------
// Datum WGS84
A  = 6378137 ;
B  = 6356752.3142 ;
e2 = 0.00669437999014 ;   
//-----------------------------------
// kleine Halbachse
B = A*Math.sqrt(1-e2) ;

// Koeffizienten
AN = (A-B)/(A+B)    ;
GP = A/(1.0 + AN)  ;
GQ = ( 1.0 + Math.pow(AN,2)/4.0 + Math.pow(AN,4)/64.0 ) ;
GR = (AN - Math.pow(AN,3)/8.0)*1.5   ;
GS = (Math.pow(AN,2)-Math.pow(AN,4)/4.0)*15.0/16.0   ;
GT = Math.pow(AN,3)*35.0/48.0  ;
GU = Math.pow(AN,4)*315.0/512.0      ;

// Bogenlänge

G = GP*(GQ*phi - GR*Math.sin(2*phi) + GS*Math.sin(4*phi) - GT*Math.sin(6*phi) + GU*Math.sin(8*phi) ) ;

return G;
}
//#######################################################################
function Ell2Gau (lamdiff,phi) {


phi     = phi * Math.PI/180 ;
lamdiff = lamdiff * Math.PI/180 ;
//----------------------------------
// Datum WGS84
A  = 6378137 ;
B  = 6356752.3142 ;
e2 = 0.00669437999014 ; 
es2=  0.00673949674228 ;
//-------------------------------------

nenner = Math.sqrt( 1.0 - e2 *( Math.pow(Math.sin(phi),2) ) ) ;
M  = A*(1.0-e2)/(Math.pow(nenner,3)) ;
N  = A/nenner ;
Rg = Math.sqrt(M*N) 
//----------------------------------------
cosB = Math.cos(phi) ;
tanB = Math.tan(phi) ;
// 2.numerische Exzentrizität
//es2= e2 /(1.-e2)  ;

eta2 = es2*Math.pow(cosB,2) ;

// Koeffizienten der Potenzreihen
a1 =   N*cosB                                             ;
a2 = - N*Math.pow(cosB,2)*tanB/2.0   ;
a3 = - N*Math.pow(cosB,3)*(1-Math.pow(tanB,2)+eta2)/ 6.0      ;
a4 =   N*Math.pow(cosB,4)*tanB*(5-Math.pow(tanB,2)+9*eta2+4*eta2*eta2)/24.0  ;
a5 =   N*Math.pow(cosB,5)*(5-18*Math.pow(tanB,2)+Math.pow(tanB,4))/120.0  ;

//--------------------------------------
// isometrische Größen
dq = 0 ;
l  = lamdiff ;
//-------------------------------------
// mathem. Bezeichnungen
dy = - a2*Math.pow(l,2) + a4* Math.pow(l,4) ;

x  =    a1*l+ a2*(2*dq*l)+a3*(3*Math.pow(dq,2)*l-Math.pow(l,3))+ a4*(4*Math.pow(dq,3)*l-4*dq*Math.pow(l,3))+ a5*(5*Math.pow(dq,4)*l-10*Math.pow(dq,2)*Math.pow(l,3)+Math.pow(l,5) ) ;
//--------------------------------------
// Meridianbogen
G = MerBogLaenge(phi) ;

//--------------------------------------
// Gauß'sche Koordinaten
y = G + dy ;
	
P = Array(x,y);
return P ;
}

//#################################################################
function  ProjectionCalc (projection,richtung,x1,y1) {

var d = 1 ;
//----------------------------
switch (richtung){
case 'for':   // geogr Koordinaten -> Projektion
    d = 1 ;
    break;
case 'inv':   // Projektion -> geogr Koordinaten 
    d = -1 ;
    break;
default:
    d = 1 ;
    break;
}
//---------------------------
x1 = x1*Math.PI/180 ;
y1 = y1*Math.PI/180 ;


switch (projection){
       
case 'Equirectangular':
    x2 = x1 ;
    y2 = y1 ;
    break;
case 'Mercator':
    if (d > 0) {
        x2 =  x1 ;
	y2 =  Math.log(Math.tan(Math.PI/4+y1/2)) ; 
    }	       
    else{
        x2 =  x1 ;
        y2  = Math.PI/2 - 2*Math.atan(Math.exp(-y1));       
    }
    break;
}
 
 
P = Array(x2,y2);
return P ;      
}