function ellipsoid_ah(a,b,fi){ //ell = ellipsoid_ah(a,,finv)
var ell=new Array();
ell.a=a;
if(fi){
ell.finv=fi;
ell.f = 1/fi;
ell.b=ell.a*(1-ell.f); // f=1-b/a => b=a(1-f)
}
if(b){
ell.b=b;
ell.f=1-ell.b/ell.a;
ell.finv = 1/ell.f;
}
ell.e2=1-(1-ell.f)*(1-ell.f);
ell.ep2=ell.e2/(1-ell.e2); // ep2=e2*a2/b2 =e2/(1-e2)
ell.es=ell.e2;
ell.e=Math.sqrt(ell.e2);
ell.ep=Math.sqrt(ell.ep2);
return ell
}
function gNormale(ell,g){
var w=Math.sqrt(1-ell.e2*Math.pow(Math.sin(g[1]),2));
var nor=ell.a/w;
return nor
}
function gRho(ell,g){
var w=Math.sqrt(1-ell.e2*Math.pow(Math.sin(g[1]),2));
return ell.a*(1-ell.e2)/w/w/w
}
function gRaz(ell,g,az){
var n=gNormale(ell,g);
var rh=gRho(ell,g);
var curb=Math.pow(Math.cos(az),2)/rh+Math.pow(Math.sin(az),2)/n;
return 1/curb
}
function enu_adh(ell,g,enu)
{
var az = Math.atan2(enu[0],enu[1]);
if (az<0){az=az+2*Math.PI}
var d = Math.sqrt(enu[0]*enu[0]+enu[1]*enu[1]);
var el = Math.atan2(enu[2],d);
var r = gRaz(ell,g,az)+g[2];
var om = Math.atan2(d,r+enu[2]);
var dell=r*om;
var dhe=(r+enu[2])/Math.cos(om)-r;
var adh=new Array();
adh=[az,dell,dhe,el,r];
return adh
}
function adh_enu(ell,g,adh)
{
var r = gRaz(ell,g,adh[0])+g[2];
var om= adh[1]/r;
var d= (r+adh[2])*Math.sin(om);
var u= (r+adh[2])*Math.cos(om)-r;
var enu=new Array();
enu[0]= d*Math.sin(adh[0]);
enu[1]= d*Math.cos(adh[0]);
enu[2]=u;
return enu
}
function dE_dN_dU(ell,g1,g2){
var dg=new Array();
for (i=0;i<3;i++){
dg[i]=g2[i]-g1[i]
}
dg[0]*=gNormale(ell,g1)*Math.cos(g1[1]);
dg[1]*=gRho(ell,g1);
return dg
}
function Inv_dE_dN_dU(ell,g1,denu){
var g2=new Array();
g2[0]=g1[0]+denu[0]*gNormale(ell,g1)*Math.cos(g1[1]);
g2[1]=g1[1]+denu[1]*gRho(ell,g1);
g2[1]=g1[2]+denu[2];
return g2
}
function geo_car(ell,g){
var x= new Array();
g[2]= g[2] ? g[2] :0;
var n=gNormale(ell,g);
x[0]=(n+g[2])*Math.cos(g[0])*Math.cos(g[1]);
x[1]=(n+g[2])*Math.sin(g[0])*Math.cos(g[1]);
x[2]=(n*(1-ell.e2)+g[2])*Math.sin(g[1]);
return x
}
/////////////////////////////////////
function car_geo(ell,x){
/////////////////////////////////////
var g= new Array();
if (x[0] == 0 && x[1] == 0) {
if (x[2]==0){
g=[0,0,-ell.a];
return g;
}
g[0]=0;
g[1]=Math.atan2(x[2],0);
g[2]=Math.abs(x[2])-ell.b;
return g
}
var p=Math.sqrt(x[0]*x[0]+x[1]*x[1]);
// var w=x[0]+p;
// g[0] = 2 * Math.atan2(x[1],w); // pas bon pour trouver 180
g[0] = Math.atan2(x[1],x[0]);
var f = Math.atan2(x[2],p*(1-ell.e2));
var h;
var i=0;
var df=999;
var fp=f;
while(df>1e-12) {
var sf=Math.sin(f);
var nor=ell.a/Math.sqrt(1-ell.e2*sf*sf);
h=p/Math.cos(f)-nor;
f=Math.atan2(x[2],p*(1-ell.e2/(1+h/nor)));
df=Math.abs(f-fp);
//document.write(f+" "+df+"
");
fp=f;
}
g[1]=f;
g[2]=h;
return g
}
//############
function rotate(ind,a,xin) {
//############
var sa=Math.sin(a);
var ca=Math.cos(a);
var xout=new Array();
if(ind==1){
xout[0]= xin[0];
xout[1]= ca*xin[1] + sa*xin[2];
xout[2]= -sa*xin[1] + ca*xin[2];
}
else if(ind==2){
xout[0]= ca*xin[0] - sa*xin[2];
xout[1]= xin[1];
xout[2]= sa*xin[0] + ca*xin[2];
}
else if(ind==3){
xout[0]= ca*xin[0]+ sa*xin[1];
xout[1]= -sa*xin[0]+ ca*xin[1];
xout[2]= xin[2];
}
return xout;
}
//##############
function gen_loc(g,xg) {
//##############
var hpi= Math.PI/2;
var xl= new Array();
xl=rotate(1,hpi-g[1],rotate(3,hpi+g[0],xg));
return xl;
}
//##############
function loc_gen(g,xl) {
//##############
var hpi= Math.PI/2;
var xg= new Array();
xg=rotate(3,-hpi-g[0],rotate(1,g[1]-hpi,xl));
return xg;
}
//#############
function tf7(t,x1) {
//#############
var x2=new Array();
var szx=x1.length;
x2[0]= x1[0]+t[0]+x1[0]*t[3] -x1[1]*t[6] +x1[2]*t[5];
x2[1]= x1[1]+t[1]+x1[1]*t[3] +x1[0]*t[6] -x1[2]*t[4];
x2[2]= x1[2]+t[2]+x1[2]*t[3] -x1[0]*t[5] +x1[1]*t[4] ;
if(szx>3){
x2[3]= x1[3]+t[7]+x1[0]*t[10] -x1[1]*t[13]+x1[2]*t[12];
x2[4]= x1[4]+t[8]+x1[1]*t[10]+x1[0]*t[13] -x1[2]*t[11];
x2[5]= x1[5]+t[9]+x1[2]*t[10]-x1[0]*t[12]+x1[1]*t[11];
}
return x2
}
/*
T08=new Array();
T08.units=[0.01, 0.01, 0.01,1e-009, 4.84813681109536e-009, 4.84813681109536e-009, 4.84813681109536e-009, 1];
T08.ITRF97=[ 0.53 , 0.01 , -4.92 , 3.37 , 0.00 , 0.00 , 0.16 , 0.01 , -0.05 , -0.32 , 0.09 , 0.00 , 0.00 , 0.02, 2005];
T08.ITRF96=[ 0.53 , 0.01 , -4.92 , 3.37 , 0.00 , 0.00 , 0.16 , 0.01 , -0.05 , -0.32 , 0.09 , 0.00 , 0.00 , 0.02, 2005];
T08.ITRF94=[ 0.53 , 0.01 , -4.92 , 3.37 , 0.00 , 0.00 , 0.16 , 0.01 , -0.05 , -0.32 , 0.09 , 0.00 , 0.00 , 0.02, 2005];
T08.ITRF93=[-3.80 , 0.19 , -5.06 , 3.86 , -2.26 , -2.43 , 0.05 , -0.28 , -0.01 , -0.24 , 0.09 , -0.11 , -0.19 , 0.07, 2005];
T08.ITRF92=[ 1.33 , 0.21 , -5.72 , 2.66 , 0.00 , 0.00 , 0.16 , 0.01 , -0.05 , -0.32 , 0.09 , 0.00 , 0.00 , 0.02, 2005];
T08.ITRF91=[ 2.53 , 1.61 , -6.32 , 4.06 , 0.00 , 0.00 , 0.16 , 0.01 , -0.05 , -0.32 , 0.09 , 0.00 , 0.00 , 0.02, 2005];
T08.ITRF90=[ 2.33 , 1.21 , -7.92 , 4.36 , 0.00 , 0.00 , 0.16 , 0.01 , -0.05 , -0.32 , 0.09 , 0.00 , 0.00 , 0.02, 2005];
T08.ITRF89=[ 2.83 , 3.61 ,-11.72 , 7.76 , 0.00 , 0.00 , 0.16 , 0.01 , -0.05 , -0.32 , 0.09 , 0.00 , 0.00 , 0.02, 2005];
T08.ITRF88=[ 2.33 , 0.01 ,-14.12 , 10.86 , 0.10 , 0.00 , 0.16 , 0.01 , -0.05 , -0.32 , 0.09 , 0.00 , 0.00 , 0.02, 2005];
T08.ITRF05=[-0.05 , -0.09 , -0.47 , 0.94 , 0.00 , 0.00 , 0.00 , 0.03 , 0.00 , 0.00 , 0.00 , 0.00 , 0.00 , 0.00, 2005];
T08.ITRF08=[ 0.00 , 0.00 , 0.00 , 0.00 , 0.00 , 0.00 , 0.00 , 0.00 , 0.00 , 0.00 , 0.00 , 0.00 , 0.00 , 0.00, 2005];
T08.ITRF00=[-0.14 , -0.12 , -1.95 , 1.74 , 0.00 , 0.00 , 0.00 , 0.01 , 0.01 , -0.18 , 0.08 , 0.00 , 0.00 , 0.00, 2005];
T08.ETRF00=[ 5.26, 4.98 ,-6.75 , 1.74 , 1.296 , 7.84 ,-12.672, 0.01 , 0.01 , -0.18 , 0.08 , 0.081 , 0.49 ,-0.792, 2005];
//T08.RGF93=[ 5.26, 4.98 ,-6.75 , 1.74 , 1.296 , 7.84 ,-12.672, 0.03 , 0.00 , -0.00 , 0.00 , 0.054 , 0.518 ,-0.782, 2005];
T08.RGF93=[ 5.26, 4.98 ,-6.75 , 1.74 , 1.296 , 7.84 ,-12.672, -0.041 , -0.022 , -0.041 , 0.00 , 0.083 , 0.534 ,-0.750, 2005];
//T08.RGF93=[ 5.30, 5.02 ,-7.47 , 2.06 , 1.62 , 9.80 ,-15.84, 0.03 , 0.00 , -0.00 , 0.00 , 0.054 , 0.518 ,-0.782, 2009];
T08.RGAF09=[ 0.07,-0.09 ,-0.47 , 0.94 , 0.00 , 0.00 , 0.00 , 0.03 , 0.00 , 0.00 , 0.00 , 0.166 , 0.650 ,-0.550, 2009];
*/
function partfo7(st_itrf,yy,dir,prefix){
//document.write(st_itrf+" "+yy+"
");
// console.log(st_itrf+" "+yy)
dir = (dir>0 || dir==undefined) ? 1 :-1;
var yyyy=yy;
// yy is THE TRUE YEAR !!
// if (yy<100){
// yyyy= (yy>=80)? 1900+yy : 2000+yy;
// }
var t=[];
var titrf=new Array();
if(prefix=='T14'){
titrf=T14[st_itrf];
}
else{
titrf=T08[st_itrf];
}
var dt=yyyy-titrf[14];
for (i=0;i<7;i++){
t[i]=(titrf[i]+titrf[i+7]*dt)*dir; //*T08.units[i]
t[i+7]=titrf[i+7]*dir; //*T08.units[i]
}
t[14]=yyyy;
return t;
}
function partfo7_u(st_itrf,yy,dir,prefix){
// console.log(st_itrf+" "+yy)
dir = dir>0 ? 1 :-1;
var yyyy=yy;
// yy is THE TRUE YEAR !!
// if (yy<100){
// yyyy= (yy>=80)? 1900+yy : 2000+yy;
// }
var t=[];
// var titrf=T08[st_itrf];
if(prefix=='T14'){
titrf=T14[st_itrf];
}
else{
titrf=T08[st_itrf];
}
var dt=yyyy-titrf[14];
for (i=0;i<7;i++){
t[i]=(titrf[i]+titrf[i+7]*dt)*dir*T14.units[i]; //
t[i+7]=titrf[i+7]*dir*T14.units[i]; //
}
t[14]=yyyy;
return t;
}
//##############
function sca_prod(x1,x2) {
//##############
var n1=x1.length;
var n2=x2.length;
var n= n1<=n2 ? n1 : n2;
var s=0;
for (i=0;i1e-20){
var sf1= tanh (lis+0.5*e*Math.log((1+e*sf)/(1-e*sf)));
dsf= Math.abs(sf1-sf);
//document.write(sf+" "+sf1+" "+dsf+"
");
sf = sf1;
}
var f= Math.asin(sf);
return f
}
//#############################################################################################################
//sub paramccf_tan{ # conformal conic parameters from def (case tangential + scale factor)
// # usage : ¶mccf_tan($ell,$l0,$f0,$k0,$x0,$y0,\%p);
//#############################################################################################################
function paramccf_tan(ell,l0,f0,k0,x0,y0){
var p= new Array();
p.ell = ell;
p.l0 = l0;
p.n = Math.sin(f0);
var r0 = gNormale(ell,[l0,f0]) * Math.cos(f0);
p.c = k0*r0/p.n * Math.exp( p.n * f_fiso(ell,f0) );
p.xs = x0;
p.ys = y0 + k0 * r0 / p.n ;
return p
}
//#############################################################################################################
//sub paramccf_sec{ # conformal conic parameters from def (case secant lat1, lat2 ...)
// # usage : ¶mccf_tan($ell,$l0,$f1,$f2,$f0,$x0,$y0,\%p);
//#############################################################################################################
function paramccf_sec(ell,l0,f1,f2,f0,x0,y0){
var p= new Array();
p.ell = ell;
p.l0 = l0;
p.f0 = f0;
var r1 = gNormale(ell,[l0,f1]) * Math.cos(f1);
var r2 = gNormale(ell,[l0,f2]) * Math.cos(f2);
p.n = ( Math.log( r2 / r1 ) ) / ( f_fiso(ell,f1) - f_fiso(ell,f2) );
p.c = r1/p.n * Math.exp( p.n * f_fiso(ell,f1) );
p.xs= x0;
p.ys= y0 + p.c * Math.exp( -p.n * f_fiso(ell,f0) );
return p
}
//#############################################################################################################
//sub confconic{ # conformal conic projection
// # usage : &confconic ($dir,$ell,\%p,\@g,\@x);
//#############################################################################################################
function confconic(dir,p,g,x,ce){
//function confconic(dir,ell,p,pt){
var gam;
var r;
if( dir> 0 ){ //# direction lon,lat -> E,N
gam = (g[0]-p.l0) * p.n;
var ll = f_fiso (p.ell,g[1]);
r = p.c * Math.exp(-(p.n*ll));
x[0] = p.xs + r * Math.sin(gam); //# E
x[1] = p.ys - r * Math.cos(gam); //# N
}
else { //# direction E,N -> lon,lat
var delx = x[0]-p.xs;
var dely = p.ys-x[1];
gam = Math.atan2(delx,dely);
g[0] = p.l0 + gam/p.n;
r = Math.sqrt( delx*delx + dely*dely );
g[1] = fiso_f(p.ell,Math.log(p.c/r)/p.n);
}
ce[0]=-gam;
ce[1]=p.n*r/gNormale(p.ell,g)/Math.cos(g[1])-1;
}
function cdir(ell){
var ca =new Array();
var ci =new Array();
var cmer =new Array();
var e2=ell.e2;
var e4=e2*e2;
var e6=e4*e2;
var e8=e6*e2;
ca[0]=1.+e2*(-1./4.+e2*( -3./64. +e2*( -5./256. +e2*( -175./16384 ))));
ca[1]= e2*(1./8. +e2*( -1./96. +e2*( -9./1024. +e2*( -901./184320 ))));
ca[2]= +e4*( 13./768. +e2*( +17./5120. +e2*( -311./737280. )));
ca[3]= +e6*( +61./15360.+e2*( +899./430080. ));
ca[4]= +e8*(+49561./41287680.);
ci[0]=ca[0];
ci[1]=e2*(1./8.+e2*(1./48. +e2*( 7./2048. +e2*( -17./184320. )))); //#DUFOUR
ci[2]= e4*(1./768. +e2*( 3./1280. +e2*( +559./368640. ))); //#NTG76 & DUFOUR
ci[3]= e6*(17./30720. +e2*( 283./430080. )); //#NTG76 & DUFOUR
ci[4]= +e8*( 4397./41287680.);
cmer[0]=ca[0];
cmer[1]= -3./8.*e2 - 3./32. *e4 -45./1024.*e6 - 105./4096. *e8;
cmer[2]= 15./256.*e4 +45./1024.*e6 + 525./16384. *e8;
cmer[3]= -35./3072.*e6 - 175./12288. *e8;
cmer[4]= 315./131072.*e8;
var tm=new Array();
tm.ca=ca;
tm.ci=ci;
tm.cmer=cmer;
return tm
}
function paramtrsv_merc(ell,l0,f0,k0,x0,y0){
p=new Array();
p.ell = ell;
p.cdir = cdir(ell);
p.l0 = l0;
p.n = k0 * ell.a;
p.xs = x0;
if (f0 != 0){
var tm=cdir(ell);
var ca=tm.ca;
var ci=tm.ci;
var cmer=tm.cmer;
var beta = cmer[0]*f0;
for (i=0;i<4;i++){
beta += cmer[i+1] * Math.sin(2*(i+1)*f0);
}
p.ys = y0 - p.n * beta;
}
else{
p.ys = y0;
}
return p
}
function utm(dir,p,g,en){
// var tm=p.cdir;
var ca=p.cdir.ca;
var ci=p.cdir.ci;
// document.write(ca+"
");
if( dir> 0 ){ //# direction lon,lat -> E,N
var dl = g[0]-p.l0;
var ll = f_fiso(p.ell,g[1]);
var sf = Math.sin(g[1]);
var w = Math.sqrt(1.- p.ell.e2*sf*sf);
var m = w/Math.cos(g[1])/cosh(ll);
var x = Math.sin(dl) / cosh(ll);
x = 0.5* Math.log((1+x)/(1-x));
m *= cosh(x);
var y = Math.atan(sinh(ll)/Math.cos(dl));
var z = [y,x];
// my zz = $ca[0]*$z;
var zz = cmul(ca[0],z);
// document.write(ca[0]+" "+z+"
"+zz+"
");
var dz = [ca[0],0];
for (i=0;i<4;i++){
zz = cadd(zz,cmul(ca[i+1],csin(cmul(2*(i+1),z))));
dz = cadd(dz,cmul(cmul(ca[i+1],ccos(cmul(2*(i+1),z))),(i+1)*2));
}
// $zz *= $$p{n};
zz = cmul(zz,p.n);
// $m *= abs($dz)*$$p{n}/$ell->a;
m *= cabs(dz)*p.n/p.ell.a;
en[0] = zz[1] + p.xs;//
en[1] = zz[0] + p.ys;
// en[2] = 0;
dz =cmul(dz, ccos(z));
//# $gam=atan2(Im($dz),Re($dz));
gam =carg(dz);
en[3]=gam/ Math.PI *180;
en[4]=(m-1)*1e6;
//# print "$$x[0] ";
}
else { //# direction E,N -> lon,lat
var xx = (en[0]-p.xs)/p.n/ca[0];
var yy = (en[1]-p.ys)/p.n/ca[0];
var zz = [yy,xx];
var z = zz;
/*
for (i=0;i<4;i++){
var zaux=cmul(-ci[i+1], csin(cmul(2*(i+1),zz)));
z = cadd(z ,zaux);
}
*/
var zi= zz;
var del=999;
var cnt=0;
while(del > 1e-20){
var zii=zz;
for (i=0;i<4;i++){
var zx=cmul(-ca[i+1]/ca[0], csin(cmul(2*(i+1),zi))); //
// $zii = $zii-$ca[$i+1]/$ca[0] * sin(2.*($i+1)*$zi);
zii = cadd(zii ,zx);
}
del=cabs(cadd(zii,cmul(-1,zi)));
zi=zii;
cnt++;
}
// document.write(cnt+" "+del+"
");
z=zi;
var x = z[1];//Im($z);
var y = z[0];//Re($z);
//# print "$z $x $y\n";
var dl = Math.atan( sinh(x)/ Math.cos(y) );
g[0] = p.l0+dl;
var ll = Math.sin(y) / cosh(x);
ll = 0.5* Math.log((1+ll)/(1-ll));
g[1] = fiso_f(p.ell,ll);
var w = Math.sqrt(1-p.ell.e2* Math.pow(Math.sin(g[1]),2));
var m = w/Math.cos(g[1])/cosh(ll);
m *= cosh(x);
// z = [Math.atan( sinh(ll)/ Math.cos(dl)),Math.sin(dl)/cosh(ll)];
var dz =[ca[0],0];
for (i=0;i<4;i++){
dz = cadd(dz,cmul(cmul(ca[i+1],ccos(cmul(2*(i+1),z))),(i+1)*2));
}
m *= (cabs(dz)*p.n/p.ell.a);
// $dz *= cos($z);
dz = cmul(dz,ccos(z));
//# $gam=atan(Im($dz)/Re($dz));
gam = carg(dz);
g[3] = gam / Math.PI *180;
g[4] = (m-1)*1e6;
}
}
function utmauto(dir,ell,g,en,zutm){
var upar=[0,0,0.9996,500000,0];
var p=new Array;
if( dir> 0 ){ //# direction lon,lat -> E,N
var zone='';
if(zutm ==''){
var lon=g[0]/Math.PI*180;
lon = lon>180? lon-360 : lon;
var fus=Math.floor((lon+180)/6)+1;
upar[0] = ( fus*6 - 183 ) / 180 * Math.PI;
var zon= fus<10 ? '0'+fus : ''+fus;
if(g[1]<0){
zone = zon.concat('S');
upar[4] = 10000000;
}
else{
zone =zon.concat('N');
}
p=paramtrsv_merc(ell,upar[0],upar[1],upar[2],upar[3],upar[4]);
}
else{
zone=zutm;
var fus= parseInt(zutm.substr(0,2),10);
var ns = zutm.substr(2,1);
upar[0] = ( fus*6 - 183 ) / 180 * Math.PI;
upar[4] = (ns=='S' || ns=='s')? 10000000 :0;
p=paramtrsv_merc(ell,upar[0],upar[1],upar[2],upar[3],upar[4]);
}
utm(1,p,g,en);
en.zone=zone;
return zone;
}
else{
// my $UZ = uc $$uz;
var fus= parseInt(zutm.substr(0,2),10);
var ns = zutm.substr(2,1);
// document.write(" zone: "+zutm+" "+fus+" "+ns+"
");
upar[0] = ( fus*6 - 183 ) / 180 * Math.PI;
upar[4] = (ns=='S' || ns=='s')? 10000000 :0;
// var p=paramtrsv_merc(ell,upar);
p=paramtrsv_merc(ell,upar[0],upar[1],upar[2],upar[3],upar[4]);
utm(-1,p,g,en);
}
}
function param_laea(ell,l0,f0,x0,y0){
var par=new Array();
par.ell=ell;
var a=par.ell.a;
var e=par.ell.e;
var e2=par.ell.e2;
var qp=(1-e2)*(1/(1-e2)-1/2/e*Math.log((1-e)/(1+e)));
var sf0=Math.sin(f0);
var q0=(1-e2)*(sf0/(1-e2*sf0*sf0)-1/2/e*Math.log((1-e*sf0)/(1+e*sf0)));
var sb0=q0/qp;
var cb0=Math.sqrt(1-sb0*sb0);
var bet0=Math.atan2(sb0,cb0); //bet0deg=bet0/$d_r;
var Rq=a*Math.sqrt(qp/2);
var m0=Math.cos(f0)/Math.sqrt(1-e2*sf0*sf0);
var D=a*m0/Rq/cb0;
par.l0=l0;
par.f0=f0;
par.x0=x0;
par.y0=y0;
par.qp=qp;
par.q0=q0;
par.beta0=bet0;
par.Rq=Rq;
par.D=D;
return par;
}
function laea(dir,par,gg,xx) {
var ell=par.ell;
var e=ell.e;
var e2=ell.e2;
var sb0=Math.sin(par.beta0);
var cb0=Math.cos(par.beta0);
if (dir > 0){ //#geo -> XY
var l=gg[0];
var f=gg[1];
var sf=Math.sin(f);
var esf=e*sf;
var q=(1-e2)*(sf/(1-esf*esf)-1/2/e*Math.log((1-esf)/(1+esf)));
var sb=q/par.qp;
var cb=Math.sqrt(1-sb*sb);
var B=par.Rq*Math.sqrt(2/(1+sb0*sb+cb0*cb*Math.cos(l-par.l0)));
xx[0]=B*par.D*cb*Math.sin(l-par.l0)+par.x0;
xx[1]=B/par.D*(cb0*sb-sb0*cb*Math.cos(l-par.l0))+par.y0;
return ;
}
else{ //#XY -> geo
var x=xx[0]-par.x0;
var y=xx[1]-par.y0;
var rho=Math.sqrt(Math.pow(x/par.D,2)+Math.pow(y*par.D,2));
if (rho==0){ gg= [par.l0,par.f0]; return }
var shce=rho/2/par.Rq;
var chce=Math.sqrt(1-shce*shce);
var ce=2*Math.atan2(shce,chce);
var sce=Math.sin(ce);
var cce=Math.cos(ce);
var q= par.qp*(cce*sb0+par.D*y*sce*cb0/rho);
if(Math.abs(q-par.qp)<1e-15){ gg= [0,Math.PI/2*q/abs(q)]; return}
gg[0]=par.l0+Math.atan2(x*sce,par.D*(rho*cb0*cce-par.D*y*sb0*sce));
var ff=Math.atan2(q/2,Math.sqrt(1-q*q/4)); //# initial lat before iterations
var eps=1e-15;
var fff;
for (i=0;i<30;i++){
var sff=Math.sin(ff);
var esff=e*sff;
fff=ff+Math.pow((1-esff*esff),2)/2/Math.cos(ff)*
(q/(1-e2)-sff/(1-esff*esff)+0.5/e*Math.log((1-esff)/(1+esff)));
if (Math.abs(fff-ff) < eps){ break}
ff=fff;
}
gg[1]=fff;
return ;
}
}
function utmauto_bis(dir,ell,g,en,zutm){
var upar=[0,0,0.9996,500000,0];
var p=new Array;
if( dir> 0 ){ //# direction lon,lat -> E,N
zone='';
var lon=g[0]/Math.PI*180;
// $lon -= 360 if $lon > 180;
lon = lon>180? lon-360 : lon;
var fus=Math.floor((lon+180)/6)+1;
// $$uz = sprintf("%02d",int(($lon+180)/6)+1);
upar[0] = ( fus*6 - 183 ) / 180 * Math.PI;
// if (fus<10){zone='0'+zone}
var zon= fus<10 ? '0'+fus : ''+fus;
// document.write(" zon : "+zon+"
");
if(g[1]<0){
zone = zon.concat('S');
upar[4] = 10000000;
}
else{
zone =zon.concat('N');
}
// zone="TM"+zone
// document.write(" zone: "+zone+"
");
// document.write(" par : "+upar+"
");
// var p=paramtrsv_merc(ell,upar);
p=paramtrsv_merc(ell,upar[0],upar[1],upar[2],upar[3],upar[4]);
utm(1,p,g,en);
en.zone=zone;
return zone;
}
else{
// my $UZ = uc $$uz;
var fus= parseInt(zutm.substr(0,2),10);
var ns = zutm.substr(2,1);
// document.write(" zone: "+zutm+" "+fus+" "+ns+"
");
upar[0] = ( fus*6 - 183 ) / 180 * Math.PI;
upar[4] = (ns=='S' || ns=='s')? 10000000 :0;
// var p=paramtrsv_merc(ell,upar);
p=paramtrsv_merc(ell,upar[0],upar[1],upar[2],upar[3],upar[4]);
utm(-1,p,g,en);
}
}
function select_ccf(ell,g){
//CCF.CC42 =paramccf_sec(GRS80, 3*DEG_RAD,(42-0.75)*DEG_RAD,(42+0.75)*DEG_RAD,42*DEG_RAD,1700000,1200000);
DEG_RAD=Math.PI/180;
RAD_DEG=180/Math.PI;
var p=new Array();
p.valid=true;
p.pn=[ell,3/180*Math.PI,0,0,0,1700000,0];
p.ps=[ell,3/180*Math.PI,0,0,0,1700000,0];
var lon=g[0]*RAD_DEG;
var lat=g[1]*RAD_DEG;
if(lat<41 || lat>51.2){
p.valid=false;
return p
}
var zn='';
var zs='';
var f0nd;
var f0sd;
if(lat<41){
f0nd=42;
f0sd=42;
}
else if(lat>51.2){
f0nd=50;
f0sd=50;
}
else if(lat>=41 && lat<43.1 && lon>8){
f0nd=42;
f0sd=42;
}
else if(lat>42 && lat<=43 && lon<8){
f0nd=43; // Math.ceil(lat);
f0sd=42; // Math.floor(lat);
zn='CC43';
zs='CC42';
}
else if(lat>50 && lat<=51.2 ){
f0nd=50;
f0sd=50;
}
else{
f0nd=Math.ceil(lat);
f0sd=Math.floor(lat);
}
f0sd= f0sd<42? f0nd : f0sd;
p.zn='CC'+f0nd;
p.zs='CC'+f0sd;
// var pdum=[ell,3/180*Math.PI,f1,f2,f0,1700000,0];
/*p.pn[2]=(f0nd-0.75)*DEG_RAD;
p.pn[3]=(f0nd+0.75)*DEG_RAD;
p.pn[4]=f0nd*DEG_RAD;
p.pn[6]=(f0nd-41)*1000000+200000;
p.ps[2]=(f0sd-0.75)*DEG_RAD;
p.ps[3]=(f0sd+0.75)*DEG_RAD;
p.ps[4]=f0sd*DEG_RAD;
p.ps[6]=(f0sd-41)*1000000+200000;
*/
p.pn=setparccf(ell,f0nd);
p.ps=setparccf(ell,f0sd);
return p
}
function select_ccf_from_coo(x){
var lat0=Math.floor(x[1]/1000000)+41;
return lat0;
}
function setparccf(ell,lat0){
var p=new Array();
p=[ell,3/180*Math.PI,0,0,0,1700000,0];
p[2]=(lat0-0.75)*Math.PI/180;
p[3]=(lat0+0.75)*Math.PI/180;
p[4]=lat0*Math.PI/180;
p[6]=(lat0-41)*1000000+200000;
return p
}
function ccfauto(dir,ell,g,xn,xs,cen,ces,zcc){
if(dir > 0){
var params=select_ccf(ell,g);
if (!params.valid){
return false
}
var pn=params.pn;
var ps=params.ps;
var parn=paramccf_sec(ell,pn[1],pn[2],pn[3],pn[4],pn[5],pn[6]);
var pars=paramccf_sec(ell,ps[1],ps[2],ps[3],ps[4],ps[5],ps[6]);
// confconic(1,ccf[1],gr,xccs,ce);
confconic(1,parn,g,xn,cen);
confconic(1,pars,g,xs,ces);
return [params.zn,params.zs]
}
else{
// var lat0=parseInt(zcc.substr(2,2),10);
var lat0=select_ccf_from_coo(xn);
var pn=setparccf(ell,lat0);
var parn=paramccf_sec(ell,pn[1],pn[2],pn[3],pn[4],pn[5],pn[6]);
confconic(-1,parn,g,xn,cen);
}
}
function DMSTodeg(dms) {
var sgn=1;
if(dms[3]=='W' || dms[3]=='w' || dms[3]=='O' || dms[3]=='o' || dms[3]=='S' || dms[3]=='s')
{sgn=-1}
return (dms[0]+dms[1]/60+dms[2]/3600)*sgn
}
function DMMTodeg(dmm) {
var sgn=1;
if(dmm[2]=='W' || dmm[2]=='w' || dmm[2]=='O' || dmm[2]=='o' || dmm[2]=='S' || dmm[2]=='s')
{sgn=-1}
return (dmm[0]+dmm[1]/60)*sgn
}
function degToDMS_arr(dec, locals,digits) {
if (digits==undefined || digits<0) {
digits= 5;
}
var positive_degrees= Math.abs(dec);
var degrees= Math.round(positive_degrees + 0.5) - 1;
var decimal_part= 60*(positive_degrees - degrees);
var minutes= Math.round(decimal_part + 0.5) - 1;
var seconds= (60*(decimal_part - minutes)).toFixed(digits);
if (seconds>=60) {
minutes= minutes+1;
seconds=(seconds-60).toFixed(digits);
}
if (minutes==60) {
degrees= degrees+1;
minutes= 0;
}
var dir= '';
dir= '' + (dec >= 0 ? locals[0] : locals[1]);
if(degrees<100 && locals[0]=='E'){
degrees="0"+degrees;
}
if(degrees<10){
degrees="0"+degrees;
}
if(minutes<10){
minutes="0"+minutes;
}
if(seconds<10){
seconds=(seconds*1).toFixed(digits);
seconds="0"+seconds;
}
var arr = [degrees,minutes,seconds,dir];
return arr;
}
function degToDMM_arr(dec, locals,digits) {
if (digits==undefined || digits<0) {
digits= 5;
}
var positive_degrees= Math.abs(dec);
var degrees= Math.round(positive_degrees + 0.5) - 1;
var minutes= (60.0*(positive_degrees - degrees)).toFixed(digits);
if (minutes==60) {
degrees= degrees+1;
minutes= 0.0.toFixed(digits);;
}
var dir= '';
dir= '' + (dec >= 0 ? locals[0] : locals[1]);
if(degrees<100 && locals[0]=='E'){
degrees="0"+degrees;
}
if(degrees<10){
degrees="0"+degrees;
}
if(minutes<10){
minutes="0"+minutes;
}
var arr = [degrees,minutes,dir];
return arr;
}
/*
//#############################################################################################################
//sub paramtrsv_mer{ # transverse Mercator parameters from def (case secant lat1, lat2 ...)
// # usage : ¶mtrsv_mer($ell,$l0,$f0,$k0,$x0,$y0,\%p);
//#############################################################################################################
my $ell = shift;
my $l0 = shift;
my $f0 = shift;
my $k0 = shift;
my $x0 = shift;
my $y0 = shift;
my $p = shift;
function paramtrsv_merc(ell,l0,f0,k0,x0,y0){
p=new Array();
p.l0 = l0;
p.n = k0 * ell.a;
p.xs = x0;
if (f0 != 0){
var tm=cdir(ell);
var ca=tm.ca;
var ci=tm.ci;
var cmer=tm.cmer;
var beta = cmer[0]*f0;
for (i=0;i<4;i++){
beta += cmer[i+1] * Math.sin(2*(i+1)*f0);
}
p.ys = y0 - p.n * beta;
}
else{
p.ys = y0;
}
return p
}
//#############################################################################################################
//sub utm{
// # usage : &utm ($dir,$ell,\@p,\@g,\@x);
//#############################################################################################################
my $dir= shift; # direction index >0 lon,lat -> E,N <=0 E,N -> lon,lat
my $ell= shift; # ref to ellipsoid object
my $p = shift; # ref to proj param array (l0, f0, k0, xs, ys)
my $g = shift; # ref to (lon,lat) rad
my $en = shift; # ref to (E,N) m
function utm(ell,p,g,en){
var tm=cdir(ell);
var ca=tm.ca;
var ci=tm.ci;
if( dir> 0 ){ //# direction lon,lat -> E,N
var dl = g[0]-p.l0;
var ll = f_fiso(ell,g[1]);
var sf=Math.sin(g[1]);
var w = sqrt(1.-ell.e2* sf*sf);
var m =w/Math.cos(g[1])/cosh(ll);
var x = Math.sin(dl) / cosh(ll);
x = 0.5* Math.log((1+x)/(1-x));
m *= cosh(x);
var y = Math.atan(sinh(ll)/Math.cos(dl));
var z = [y,x];
// my zz = $ca[0]*$z;
var zz = cmul(ca[0],z];
var dz = [ca[0],0];
for (i=0;i<4;i++){
zz = cadd(zz,cmul(ca[i+1],csin(cmul(2*(i+1),z))));
dz = cadd(dz,cmul(cmul(ca[i+1],ccos(cmul(2*(i+1),z))),(i+1)*2));
}
// $zz *= $$p{n};
zz = cmul(zz,p.n};
// $m *= abs($dz)*$$p{n}/$ell->a;
m = cabs($dz)*p.n/ell.a;
en[0] = zz[1] + p.xs;
en[1] = zz[0] + p.ys;
dz =cmul(dz, ccos(z));
//# $gam=atan2(Im($dz),Re($dz));
gam =carg(dz);
//# print "$$x[0] ";
}
else { //# direction E,N -> lon,lat
var xx = (en[0]-p.xs)/p.n/ca[0];
var yy = (en[1]-p.ys)/p.n/ca[0];
var zz = [yy,xx];
var z = zz;
for (my i=0;i<4;i++){
var zaux=cmul(-ci[i+1], csin(cmul(2*(i+1),zz));
z = cadd(z ,zaux);
}
var zi= zz;
//# for(my $k=0;$k<15;$k++){
var del=1;
while(del>1e-12){
var zii=zz;
for (my i=0;i<4;i++){
var zaux=cmul(-ca[i+1]/ca[0], csin(cmul(2*(i+1),zi));
// $zii = $zii-$ca[$i+1]/$ca[0] * sin(2.*($i+1)*$zi);
zii = cadd(zii ,zaux);
}
var del=cabs(zii-zi);
zi=zii;
}
z=zi;
var x = z[1];//Im($z);
var y = z[0];//Re($z);
//# print "$z $x $y\n";
var dl = atan( sinh(x)/ Math.cos(y) );
g[0] = p.l0+dl;
var ll = Math.sin(y) / cosh(x);
ll = 0.5* Math.log((1+ll)/(1-ll));
g[1] = fiso_f(ell,ll);
var w = Math.sqrt(1-ell.e2* Math.pow(Math.sin(g[1]),2));
var m = w/Math.cos(g[1])/cosh(ll);
m *= cosh(x);
z = [Math.atan( sinh(ll)/ Math.cos(dl)),Math.sin(dl)/cosh(ll)];
var dz =[ca[0],0];
for (i=0;i<4;i++){
dz = cadd(dz,cmul(cmul(ca[i+1],ccos(cmul(2*(i+1),z))),(i+1)*2));
}
m *= (cabs(dz)*p.n/ell.a);
// $dz *= cos($z);
dz = cmul(dz,ccos($z));
# $gam=atan(Im($dz)/Re($dz));
gam = carg(dz);
g[3] = gam / Math.PI *180;
g[4] = (m-1)*1e6;
}
}
//sub cdir{
function cdir(ell){
var ca =new Array();
var ci =new Array();
var cmer =new Array();
var e2=ell.e2;
var e4=e2*e2;
var e6=e4*e2;
var e8=e6*e2;
ca[0]=1.+e2*(-1./4.+e2*( -3./64. +e2*( -5./256. +e2*( -175./16384 ))));
ca[1]= e2*(1./8. +e2*( -1./96. +e2*( -9./1024. +e2*( -901./184320 ))));
ca[2]= +e4*( 13./768. +e2*( +17./5120. +e2*( -311./737280. )));
ca[3]= +e6*( +61./15360.+e2*( +899./430080. ));
ca[4]= +e8*(+49561./41287680.);
ci[0]=ca[0];
ci[1]=e2*(1./8.+e2*(1./48. +e2*( 7./2048. +e2*( -17./184320. )))); //#DUFOUR
ci[2]= e4*(1./768. +e2*( 3./1280. +e2*( +559./368640. ))); //#NTG76 & DUFOUR
ci[3]= e6*(17./30720. +e2*( 283./430080. )); //#NTG76 & DUFOUR
ci[4]= +e8*( 4397./41287680.);
cmer[0]=ca[0];
cmer[1]= -3./8.*e2 - 3./32. *e4 -45./1024.*e6 - 105./4096. *e8;
cmer[2]= 15./256.*e4 +45./1024.*e6 + 525./16384. *e8;
cmer[3]= -35./3072.*e6 - 175./12288. *e8;
cmer[4]= 315./131072.*e8;
var tm=new Array();
tm.ca=ca;
tm.ci=ci;
tm.cmer=cmer;
return tm
}
#############################################################################################################
sub utmauto{
# usage : &utmauto ($dir,$ell,\$uz,\@g,\@en);
#############################################################################################################
my $dir = shift; # direction index >0 lon,lat -> E,N <=0 E,N -> lon,lat
my $ell = shift; # ref to ellipsoid object
my $uz = shift ; # ref to UTM zone "ffX" (ff:zone (01..60) X:N or S) in for E,N -> lon,lat out for lon,lat -> E,N
my $g = shift; # ref to (lon,lat) rad
my $en = shift; # ref to (E,N) m
my $pi = 4* atan(1);
function utmauto(dir,ell,zone,g,en){
var utmpar=[0,0,0.9996,500000,0];
if( dir> 0 ){ //# direction lon,lat -> E,N
var lon=g[0]/Math.PI*180;
// $lon -= 360 if $lon > 180;
lon = lon>180? lon-360 ; lon;
var zone=Math.floor((lon+180)/6);
// $$uz = sprintf("%02d",int(($lon+180)/6)+1);
utmpar[0] = ( zone*6 - 183 ) / 180 * Math.PI;
if (zone<10){zone='0'+zone}
if(g[1]<0){
zone += 'S';
utmpar[4] = 10000000;
}
else{
zone += 'N';
}
var p=paramtrsv_mer(ell,utmpar);
utm(1,ell,p,g,en);
}
else{
// my $UZ = uc $$uz;
fus= parseInt(substr(zone,0,2),10);
var ns = substr(zone,2,1);
utmpar[0] = ( fus*6 - 183 ) / 180 * Math.PI;
utmpar[4] = ns=='S'? 10000000 :0;
var p=paramtrsv_mer(ell,@utmpar);
utm(-1,ell,p,g,en);
}
}
sub param_laea{
my $ell=shift;
my $l0=shift;
my $f0=shift;
my $x0=shift;
my $y0=shift;
function param_laea(ell,l0,f0,x0,y0)
var par=new Array();
var a=ell.a;
var e=ell.e;
var e2=ell.e2;
var qp=(1-e2)*(1/(1-e2)-1/2/e*Math.log((1-e)/(1+e)));
var sf0=Math.sin(f0);
var q0=(1-e2)*(sf0/(1-e2*sf0*sf0)-1/2/e*Math.log((1-e*sf0)/(1+e*sf0)));
var sb0=q0/qp;
var cb0=Math.sqrt(1-sf0*sf0);
var bet0=Math.atan2(sb0,cb0); //bet0deg=bet0/$d_r;
var Rq=a*Math.sqrt(qp/2);
var m0=Math.cos(f0)/Math.sqrt(1-e2*sf0*sf0);
var D=a*m0/Rq/cb0;
par.ell=ell;
par.l0=l0;
par.f0=f0;
par.x0=x0;
par.y0=y0;
par.qp=qp;
par.q0=q0;
par.beta0=bet0;
par.Rq=Rq;
par.D=D;
return par;
}
sub laea{
my $par=shift; #ref to hash of parameters
my $gg=shift; #ref to coord. vector ( [lon,lat] )
my $xx=shift; #ref to coord. vector ( [x,y] )
my $dir=shift; #transformation direction (>0 : [lon,lat] -> [x,y] ; inverse othewise)
function laea(dir,par,gg,xx) {
var e=par.ell.e;
var e2=par.ell.e2;
var sb0=Math.sin(par.beta0);
var cb0=Math.cos(par.beta0);
if (dir > 0){ //#geo -> XY
var l=gg[0];
var f=gg[1];
var sf=Math.sin(f);
var esf=e*sf;
var q=(1-e2)*(sf/(1-esf*esf)-1/2/e*Math.log((1-esf)/(1+esf)));
var sb=q/par.qp;
var cb=sqrt(1-sb*sb);
var B=par.Rq*Math.sqrt(2/(1+sb0*sb+cb0*cb*Math.cos(l-par.l0)));
xx[0]=B*par.D*cb*Math.sin(l-par.l0)+par.x0;
xx[1]=B/par.D*(cb0*sb-sb0*cb*cos(l-par.l0))+par.y0;
return ;
}
else{ //#XY -> geo
var x=xx[0]-par.x0;
var y=xx[1]-par.y0;
var rho=Math.sqrt(Math.pow(x/par.D,2)+Math.pow(y*par.D,2));
if (rho==0){ gg= [par.l0,par.f0]; return }
var shce=rho/2/par.Rq;
var chce=Math.sqrt(1-shce*shce);
var ce=2*Math.atan2(shce,chce);
var sce=Math.sin(ce);
var cce=Math.cos(ce);
var q= rho? par.qp*(cce*sb0+par.D*y*sce*cb0/rho): ;
if(abs(q-par.qp)<1e-15){ gg= [0,Math.PI/2*q/abs(q)]; return}
gg[0]=par.l0+Math.atan2(x*sce,par.D*(rho*cb0*cce-par.D*y*sb0*sce));
var ff=Math.atan2(q/2,sqrt(1-q*q/4)); //# initial lat before iterations
var eps=1e-15;
var fff;
for (i=0;i<30;i++){
var sff=Math.sin(ff);
var esff=e*sff;
fff=ff+Math.pow((1-esff*esff),2)/2/Math.cos(ff)*
(q/(1-e2)-sff/(1-esff*esff)+0.5/e*Math.log((1-esff)/(1+esff)));
//#print "$i $ff $fff\n";
if (abs(fff-ff) < eps){ break}
ff=fff;
}
gg[1]=fff;
return ;
}
}
*/