Back to: Main Page > Mathematical programs

Program of cartographical mappings
Sources of the mathematical module

Here are the functions used to compute the diverse projections.

Some notations:

The first two functions (AzTrans and CylTrans) are the coordinates changes to move the system of coordinates of the sphere
lo=longitude, la=latitude
PProj is a structure grouping all the necessary information
For azimuthal and conical projections, k is the generalised Aitoff transformation ratio (divide longitudes by k, apply the projections, multiply the abscissae by k) and miller is the same thing for latitudes and ordinates. For cylindrical projections, k is an affinity ratio (to get some automecoic parallel) computed in the function SetParams at the end of the file.
For conical projections, laco1 and laco2 are the latitudes of the two automecoical parallels ; y2, cco and nco are parameters obtained from the function SetParams at the end of the file.
In the SetParams function, pclo, pcla, ppc2lo, pc2la, pc3lo, pc3la are the longitude and latitude of the first, second and third point of construction respectively. Generally, the first point specifies the origin of the coordinate system on the sphere; the second point specifies the direction of the system of coordinates, and optionnally an automecoical parallel; the third point may specify a second automecoical parallel.

I'm sorry for the typesetting of the code (especially indentation), which has suffered from several electronic translations and my naturally obfuscated writing style...


float sqr(float a){return a*a;}

void AzTrans(PTWProj PProj,float lo,float la,float*x,float*y)
{if(PProj->pcla==pi/2){
	*x=lo-PProj->az;if(*x>pi)*x-=2*pi;else if(*x<-pi)*x+=2*pi;
	*y=la;
	return;};
*x=lo-PProj->pclo;
float cp=cos(*x),ct=cos(la),c0=cos(PProj->pcla),sp=sin(*x),st=sin(la),s0=sin(PProj->pcla);
if(*x!=0||la!=PProj->pcla)*x=atan2(sp*ct,cp*ct*s0-st*c0);else *x=0;
*y=cp*ct*c0+st*s0;if(fabs(*y)>1)*y=sig(*y);
*y=asin(*y);
*x-=PProj->az;if(*x>pi)*x-=2*pi;else if(*x<-pi)*x+=2*pi;
}

void CylTrans(PTWProj PProj,float lo,float la,float*x,float*y)
{if(PProj->pcla==0.0&&PProj->az==0.0){
	*x=lo-PProj->pclo;if(*x>pi)*x-=2*pi;else if(*x<-pi)*x+=2*pi;
	*y=la;
	return;};
int j1,j2;
float m1[3],m2[3],lc2=cos(la);
m1[0]=lc2*cos(lo);m1[1]=lc2*sin(lo);m1[2]=sin(la);
for(j1=0;j1<3;j1++){
	m2[j1]=0;
	for(j2=0;j2<3;j2++)m2[j1]+=PProj->ma[j1][j2]*m1[j2];
  };
if(fabs(m2[2])>1)m2[2]=sig(m2[2]);
*y=asin(m2[2]);
*x=atan2(m2[1],m2[0]);
}

void Sanson(PTWProj PProj,float lo,float la,float*x,float*y)
{*x=lo*cos(la/PProj->miller);
*y=la*PProj->k;
}

void Mercator(PTWProj PProj,float lo,float la,float*x,float*y)
{float l=la/PProj->miller/2+pi/4;
if(fabs(l)==pi/2||tan(l)<=0){*x=err;return;};
*x=lo;
*y=log(tan(l))*PProj->k;
}

void CylLambert(PTWProj PProj,float lo,float la,float*x,float*y)
{*x=lo;
*y=sin(la/PProj->miller)*PProj->k;
}

void CylCentrale(PTWProj PProj,float lo,float la,float*x,float*y)
{if(fabs(la/=PProj->miller)==pi/2){*x=err;return;};
*x=lo;
*y=tan(la)*PProj->miller;
}

void Gnomonique(PTWProj PProj,float lo,float la,float*x,float*y)
{if(fabs(lo/=PProj->aitoff)>=pi/2||fabs(la/=PProj->miller)==pi/2){*x=err;return;};
*x=tan(lo)*PProj->aitoff;
*y=tan(la)/cos(lo)*PProj->miller;
}

void Orthographique(PTWProj PProj,float lo,float la,float*x,float*y)
{if(fabs(lo/=PProj->aitoff)>pi/2){*x=err;return;};
*x=sin(lo)*cos(la/=PProj->miller)*PProj->aitoff;
*y=sin(la)*PProj->miller;
}

void Stereographique(PTWProj PProj,float lo,float la,float*x,float*y)
{float a=cos(lo/=PProj->aitoff),b=cos(la/=PProj->miller);if(a*b==-1){*x=err;return;};
*x=2*sin(lo)*b/(1+a*b)*PProj->aitoff;
*y=2*sin(la)/(1+a*b)*PProj->miller;
}

void CylGall(PTWProj PProj,float lo,float la,float*x,float*y)
{*x=lo;
*y=(1+M_SQRT2)*tan(la/PProj->miller/2)*PProj->k;
}

void AzLambert(PTWProj PProj,float lo,float la,float*x,float*y)
{float d=sqrt(2/(1+cos(lo/=PProj->aitoff)*cos(la/=PProj->miller)));
*x=d*sin(lo)*cos(la)*PProj->aitoff;
*y=d*sin(la)*PProj->miller;
}

void AzHatt(PTWProj PProj,float lo,float la,float*x,float*y)
{if(lo==0&&la==0){*x=0;*y=0;return;};
float l1,l2,ct=cos(la/=PProj->miller);
l1=atan2(sin(lo/=PProj->aitoff)*ct,-sin(la));
l2=acos(cos(lo)*ct);
*x=l2*sin(l1)*PProj->aitoff;
*y=-l2*cos(l1)*PProj->miller;
}

void Mollweide(PTWProj PProj,float lo,float la,float*x,float*y)
{int n;float k=la/=PProj->miller,l=pi*sin(la);
if(fabs(la)!=pi/2)for(n=1;n<5;n++)k=(2*k*cos(2*k)-sin(2*k)+l)/(2+2*cos(2*k));
*x=2*sqrt(2)/pi*lo*cos(k);
*y=sin(k)*sqrt(2)*PProj->miller;
}

void ConLambert(PTWProj PProj,float lo,float la,float*x,float*y)
{if(la==-pi/2){*x=err;return;};
float a;if(la==pi/2)a=0;else a=exp(log(tan((pi/2-la)/2))*PProj->nco)/PProj->k;
*x=a*sin(lo*=PProj->nco);*y=-a*cos(lo);
}

void ConCentrale(PTWProj PProj,float lo,float la,float*x,float*y)
{float a;
if(fabs(la-PProj->cco)>=pi/2){*x=err;return;};
a=PProj->k*tan(la-PProj->cco)+PProj->y2;
*x=-a*sin(lo*=PProj->nco);*y=a*cos(lo);
}

void ConEquiv(PTWProj PProj,float lo,float la,float*x,float*y)
{float a=sqrt(PProj->cco-2/PProj->nco*sin(la));
*x=a*sin(lo*=PProj->nco);*y=-a*cos(lo);
}

void Bonne(PTWProj PProj,float lo,float la,float*x,float*y)
{float a=PProj->y2+(PProj->laco1-la),nco=cos(la)/a;
*x=a*sin(lo*=nco);*y=-a*cos(lo);
}

void Delisle(PTWProj PProj,float lo,float la,float*x,float*y)
{float a=PProj->y2+PProj->laco1-la;
*x=a*sin(lo*=PProj->nco);*y=-a*cos(lo);
}

void Polyconique(PTWProj,float lo,float la,float*x,float*y)
{if(la!=0.0){float s=sin(la),c=cos(la);
	*x=c/s*sin(s*lo);
	*y=la-pi/2-c/s*(cos(s*lo)-1);
  }else{
	*x=lo;*y=-pi/2;
  };
}

void Equid2pts(PTWProj PProj,float lo,float la,float*x,float*y)
{float a=acos(cos(lo)*cos(la)),b=acos(cos(lo-PProj->k)*cos(la));
*x=(a*a-b*b)/(2*PProj->k)+PProj->k/2;
*y=sqrt(fabs(a*a-*x*(*x)))*sig(la);
}

void Azimutale(PTWProj PProj,float lo,float la,float*x,float*y)
{float v=PProj->v,a=cos(lo/=PProj->aitoff),b=cos(la/=PProj->miller);
if(a*b<=-v||(v!=0.0&&a*b<-1/v)){*x=err;return;};
*x=(v+1)*sin(lo)*b/(v+a*b)*PProj->aitoff;
*y=(v+1)*sin(la)/(v+a*b)*PProj->miller;
}

void Cylindrique(PTWProj PProj,float lo,float la,float*x,float*y)
{float v=PProj->v,a=cos(la/=PProj->miller)+v;
if(a<=0.0){*x=err;return;};
*x=lo;*y=(1+v)/a*sin(la)*PProj->k;
}

void Breusing(PTWProj PProj,float lo,float la,float*x,float*y)
{if(lo==0&&la==0){*x=0;*y=0;return;};
float l1,d,ct=cos(la/=PProj->miller),r;
l1=atan2(sin(lo/=PProj->aitoff)*ct,-sin(la));
r=acos(cos(lo)*ct);
d=2*sqrt(sin(r/2)*tan(r/2));
*x=d*sin(l1)*PProj->aitoff;
*y=-d*cos(l1)*PProj->miller;
}

void Breusing2(PTWProj PProj,float lo,float la,float*x,float*y)
{if(lo==0&&la==0){*x=0;*y=0;return;};
float l1,d,ct=cos(la/=PProj->miller),r;
l1=atan2(sin(lo/=PProj->aitoff)*ct,-sin(la));
r=acos(cos(lo)*ct);
d=4*tan(r/4);
*x=d*sin(l1)*PProj->aitoff;
*y=-d*cos(l1)*PProj->miller;
}

void Littrow(PTWProj PProj,float lo,float la,float*x,float*y)
{if(fabs(la/=PProj->miller)==pi/2||fabs(lo/=PProj->aitoff)>pi/2){*x=err;*y=err;return;};
*x=sin(lo)/cos(la)*PProj->aitoff;
*y=cos(lo)*tan(la)*PProj->miller;
}

void Apianus(PTWProj PProj,float lo,float la,float*x,float*y)
{*x=lo/pi*sqrt(pi*pi-4*sqr(la/PProj->miller));
*y=la;
}

void Ortelius(PTWProj PProj,float lo,float la,float*x,float*y)
{*x=lo/2*(1+sqrt(pi*pi-4*sqr(la/PProj->miller))/pi);
*y=la;
}

void Winkel(PTWProj PProj,float lo,float la,float*x,float*y)
{float d=sqrt(2/(1+cos(lo/=PProj->aitoff)*cos(la/=PProj->miller)));
*x=(d*sin(lo)*cos(la)+lo)/2*PProj->aitoff;
*y=(d*sin(la)+la*PProj->k)/2*PProj->miller;
}

void Eckert4(PTWProj PProj,float lo,float la,float*x,float*y)
{float k=(la/=PProj->miller),r=2/sqrt(pi*(4+pi)),l=(4+pi)*sin(la);
if(fabs(k)!=pi/2)for(int n=1;n<5;n++)
	k-=(2*k+4*sin(k)+sin(2*k)-l)/(2+4*cos(k)+2*cos(2*k));
*x=r*(1+cos(k))*lo;
*y=r*pi*sin(k)*PProj->miller;
}

void Collignon(PTWProj PProj,float lo,float la,float*x,float*y)
{*y=sqrt(pi)*(1-sqrt(1-fabs(sin(la/PProj->miller))))*PProj->miller*sig(la);
*x=2/pi*lo*(sqrt(pi)-fabs(*y/PProj->miller));
}

void Braun(PTWProj PProj,float lo,float la,float*x,float*y)
{*x=lo;
*y=7*sin(la/=PProj->miller)/(2+5*cos(la))*PProj->k;
}

void Globulaire(PTWProj PProj,float lo,float la,float*x,float*y)
{if((lo/=PProj->aitoff)==0.0){*x=0;*y=la;return;};
if(fabs(la/=PProj->miller)==pi/2){*x=0;*y=pi/2*PProj->miller;return;};
if(la==0.0){*x=lo*PProj->aitoff;*y=0;return;};
float s=sig(la);la*=s;
float cla=pi/2-la,d,e,g,r,yc,m2,r2,xc,l,al,th;
d=sqrt(la*la+pi*pi/4-la*pi*cos(cla));
e=pi-2*asin(pi/2*sin(cla)/d);
g=e+cla;
r=pi/2*sin(cla)/sin(e);
yc=pi/2*sin(g)/sin(e);
m2=2*atan(2*lo/pi);
r2=pi/2/sin(m2);xc=pi/2*cos(m2)/sin(m2);
l=asin(cos(g)*sin(m2));
al=atan2(xc,yc);
th=asin(sin(al)*cos(l)/cos(m2))-al;
*x=r*sin(th)*PProj->aitoff;
*y=r2*sin(th+l)*PProj->miller*s;
}

void Polyedrique(PTWProj PProj,float lo,float la,float*x,float*y)
{if(fabs(lo/=PProj->aitoff)>=pi/2){*x=err;*y=err;return;};
*x=(1-sin(PProj->k)*sin(la/=PProj->miller))/cos(PProj->k)*sin(lo)/cos(lo)*PProj->aitoff;
*y=(sin(la)-sin(PProj->k))/cos(PProj->k)*PProj->miller*cos(PProj->k)/cos(PProj->k*PProj->miller);
}

void Trapeziforme(PTWProj PProj,float lo,float la,float*x,float*y)
{*y=la;
*x=lo*(PProj->nco*la+PProj->cco);
}

void BalanceErr(PTWProj PProj,float lo,float la,float*x,float*y)
{if(lo==0&&la==0){*x=0;*y=0;return;};
float l1,l2,d,ct=cos(la/=PProj->miller);
l1=atan2(sin(lo/=PProj->aitoff)*ct,-sin(la));
l2=asin(cos(lo)*ct);
l2=pi/2-l2;if((l2/=2)==pi/2){*x=err;*y=err;return;};
if(l2==0){*x=0;*y=0;return;};
d=(-2/tan(l2)*log(cos(l2))+tan(l2)*PProj->k)*2/(PProj->k+1);
*x=d*sin(l1)*PProj->aitoff;
*y=-d*cos(l1)*PProj->miller;
}

void Lagrange(PTWProj PProj,float lo,float la,float*x,float*y)
{if(fabs(lo/=PProj->aitoff)*PProj->v>=pi/2){*x=err;*y=err;return;};
if((la/=PProj->miller)==pi/2){*x=0;*y=PProj->miller/PProj->v;return;};
if(la==-pi/2){*x=0;*y=-PProj->miller/PProj->v;return;};
float u=log(tan(pi/4-la/2));complex u1=complex(0,-1)*(u+PProj->k);
complex z=tan(PProj->v*(u1+lo))/PProj->v;
*x=real(z)*PProj->aitoff;
*y=imag(z)*PProj->miller;
}

void PersoProj(PTWProj PProj,float lo,float la,float*x,float*y)
{*PProj->lovar->pval=(lo/=PProj->aitoff);
*PProj->lavar->pval=(la/=PProj->miller);
if(PProj->zPersoProj){
	float l=la/2+pi/4;
	if(l==0||l==pi/2){*x=err;*y=err;return;};
  *PProj->zvar->pval=complex(lo,log(tan(l)));
	};
complex p=PProj->PersoProjExpr.Val();
if(p==complex(ErrVal)){*x=err;*y=err;return;};
*x=real(p)*PProj->aitoff;*y=imag(p)*PProj->miller;
}

BOOL IsCyl(void(*projeq)(PTWProj,float,float,float*,float*))
{return(projeq==Quadratique||projeq==CylCentrale||projeq==Mercator||
	projeq==CylLambert||projeq==Sanson||projeq==Cylindrique||
	projeq==Mollweide||projeq==CylGall||projeq==Equid2pts||
	projeq==Eckert4||projeq==Braun);}

BOOL IsAz(void(*projeq)(PTWProj,float,float,float*,float*))
{return(projeq==Stereographique||projeq==AzHatt||projeq==AzLambert
	||projeq==Gnomonique||projeq==Orthographique||projeq==Azimutale
	||projeq==Breusing||projeq==Breusing2||projeq==Littrow
	||projeq==Apianus||projeq==Winkel||projeq==Ortelius
	||projeq==Collignon||projeq==Globulaire||projeq==Polyedrique
	||projeq==BalanceErr||projeq==Lagrange);}

BOOL IsCon(void(*projeq)(PTWProj,float,float,float*,float*))
{return(projeq==ConLambert||projeq==ConCentrale||projeq==ConEquiv
	||projeq==Bonne||projeq==Delisle||projeq==Polyconique
	||projeq==Trapeziforme);}

void TWProj::SetParams()
{if(projeq==NULL)return;
float x,y,laco2;
if(projeq==Equid2pts&&pcla==pc2la&&pclo==pc2lo)pcla+=0.01;
if(IsAz(projeq)){
        az=0;if(projeq==BalanceErr)k=0;
        if(pclo!=pc2lo||pcla!=pc2la){
                AzTrans(this,pc2lo,pc2la,&x,&y);
                az=-x;if(projeq==Polyedrique)az+=pi;
                if(projeq==BalanceErr){k=pi/2-y;k=1/sqr(tan(k/2))*log(1/sqr(cos(
k/2)));};
                };
        if(aitoff==0.0)aitoff=1;if(miller==0.0)miller=1;
        };
if(IsCyl(projeq)||(projeq==PersoProj&&!PersoProjTrans)){
        az=0;
        if(pclo!=pc2lo||pcla!=pc2la){
                AzTrans(this,pc2lo,pc2la,&x,&y);
          az=pi/2-x;
                };
        };
if(IsCon(projeq)||(projeq==PersoProj&&PersoProjTrans)){
        az=0;
        if(pclo!=pc2lo||pcla!=pc2la){
                AzTrans(this,pc2lo,pc2la,&x,&y);
                az=x;if(projeq==Trapeziforme)az=-az-pi;
                };
        AzTrans(this,pc2lo,pc2la,&x,&y);laco1=y;
        AzTrans(this,pc3lo,pc3la,&x,&y);laco2=y;
        if(projeq==Polyconique)return;
        if((laco1==laco2&&projeq!=ConCentrale)||projeq==Bonne){
                nco=sin(laco1);
                if(projeq==Bonne||projeq==Delisle){
                        if(laco1==0.0)laco1=0.0001;if(nco==0)nco=0.0001;
                        y2=cos(laco1)/sin(laco1);
      };
                if(projeq==ConEquiv){
                        if(laco1==0.0)laco1=0.0001;if(nco==0)nco=0.0001;
                        cco=cos(laco1)/sin(laco1);cco=2+cco*cco;
                        };
                if(projeq==ConLambert){
                        y2=0;k=1;projeq(this,0,laco1,&x,&y2);
                        k=-nco*y2/cos(laco1);if(k==0.0)k=1;
                        };
                }else{
                if(projeq==ConCentrale){
                        if(!((laco1+laco2)/2))laco1+=0.0001;
      cco=(laco1+laco2)/2;
                        nco=sin(cco);
                        k=cos((laco1-laco2)/2);
                        y2=-k/tan(cco);
                        };
                if(projeq==ConEquiv){
            nco=(sqr(cos(laco1))-sqr(cos(laco2)))/2/(sin(laco2)-sin(laco1));
          cco=(sqr(cos(laco1))+2*nco*sin(laco1))/sqr(nco);
                        };
                if(projeq==Delisle){
            y2=(laco2-laco1)/(1-cos(laco2)/cos(laco1));
                        nco=cos(laco1)/y2;
            };
                if(projeq==ConLambert){
                        nco=(log(sin(pi/2-laco1))-log(sin(pi/2-laco2)))/(log(tan
((pi/2-laco1)/2))-log(tan((pi/2-laco2)/2)));
                        y2=0;k=1;projeq(this,0,laco1,&x,&y2);
                        k=-nco*y2/cos(laco1);if(k==0.0){y2=0;k=1;projeq(this,0,l
aco2,&x,&y2);k=-nco*y2/cos(laco2);};
                        };
                };
        };
if(trans==CylTrans){
        cp=cos(pclo);sp=sin(pclo);
        ct=cos(pcla);st=sin(pcla);
        ca=cos(az);sa=sin(az);
        ma[0][0]=ct*cp;ma[0][1]=sp*ct;ma[0][2]=st;
        ma[1][0]=cp*st*sa-sp*ca;
        ma[1][1]=cp*ca+sp*st*sa;ma[1][2]=-ct*sa;
        ma[2][0]=-cp*st*ca-sp*sa;
        ma[2][1]=cp*sa-sp*st*ca;ma[2][2]=ct*ca;
        CylTrans(this,pc3lo,pc3la,&x,&y);
        if(cos(y)==0.0)y=0;
        if(projeq==Quadratique)k=1/cos(y);
        if(projeq==Mercator)k=cos(y/miller)/cos(y)*miller;
        if(projeq==CylGall)k=2*sqr(cos(y/miller/2))/(1+sqrtl(2))/cos(y)*miller;
        if(projeq==CylCentrale)k=sqr(cos(y/miller))/cos(y)*miller;
        if(projeq==CylLambert)k=1/cos(y/miller)/cos(y)*miller;
        if(projeq==Equid2pts){CylTrans(this,pc2lo,pc2la,&x,&y);k=x;};
        if(projeq==Cylindrique)k=1/cos(y)*sqr(v+cos(y/miller))/(1+v)/(1+v*cos(y/
miller))*miller;
        if(projeq==Sanson)k=cos(y/miller)/cos(y);
        if(projeq==Winkel)k=1/cos(y);
        if(projeq==Braun)k=sqr(5*cos(y/miller)+2)/(14*cos(y/miller)+35)/cos(y)*m
iller;
        if(projeq==Trapeziforme){
                CylTrans(this,pc2lo,pc2la,&x,&y);laco1=y;
                CylTrans(this,pc3lo,pc3la,&x,&y);laco2=y;
                if(laco1!=laco2){nco=(cos(laco2)-cos(laco1))/(laco2-laco1);cco=c
os(laco1)-nco*laco1;}
                        else{nco=-sin(laco1);cco=cos(laco1)-nco*laco1;};
                };
        if(projeq==Polyedrique){
                CylTrans(this,pc2lo,pc2la,&x,&y);k=y/miller;
                };
        if(projeq==Lagrange){
                CylTrans(this,pc2lo,pc2la,&x,&y);if(fabsl(y/2)==pi/4)y*=0.9999;
                k=-log(tan(pi/4-y/2));
                CylTrans(this,pc3lo,pc3la,&x,&y);v=sqrt(1+sqr(cos(y)))/2;
        };
        };
}

Back to: Main Page > Mathematical programs

To leave a comment: contact (domain) yann-ollivier.org

Program of cartographical mappingsSources of the mathematical module

Program of cartographical mappings
Sources of the mathematical module