Skip to content

GitLab

Dear CS10-Gitlab-users, on Thursday, Feb 3 there will be maintenance. That will lead to a downtime of the CS10-Gitlab-service including Subversion and Mattermost chat from 09:30. This might take the whole day since we don't know how long it is going to take. We are sorry for the inconvenience! Best regards, CS10-Admin-Team

Commit 26f0d171 authored by Phillip Lino Rall's avatar Phillip Lino Rall
Browse files

increased efficiency of iterpolator. trilinear interpolation not worth it, as... 

increased efficiency of iterpolator. trilinear interpolation not worth it, as it increases computational effort quiet a lot, while barely being more accurate
parent 01bc6394
Hide whitespace changes
Inline Side-by-side
Showing with 492 additions and 215 deletions
program/source/interpol/interpol.cc
View file @ 26f0d171
......@@ -53,10 +53,11 @@ bool contained_in_tet(D3vector lam) {
}
bool contained_in_tet_strong(D3vector lam) {
if(lam.x < -0.2) return false;
if(lam.y < -0.2) return false;
if(lam.z < -0.2) return false;
if(lam.x + lam.y + lam.z > 1.2) return false;
double limit = 0.2; // 0.2
if(lam.x < -limit) return false;
if(lam.y < -limit) return false;
if(lam.z < -limit) return false;
if(lam.x + lam.y + lam.z > 1+limit) return false;
return true;
}
......@@ -65,8 +66,10 @@ bool new_lam_better(D3vector lam_old, D3vector lam_new) {
if (MIN(lam_new ) < -0.2 || MAX(lam_new) > 1.2) return false;
if (MIN(lam_old ) < -0.2 || MAX(lam_old) > 1.2) return true;
if (MIN(lam_new) < MIN(lam_old)) return true;
if( (MIN(lam_new ) < 0 || MAX(lam_new) >1 )&& (MIN(lam_old ) >= 0 || MAX(lam_old) <=1 )) return false;
if (MIN(lam_new) > MIN(lam_old)) return true;
if (MAX(lam_new) < MAX(lam_old)) return true;
return false;
}
bool new_lam_worse(D3vector lam_old, D3vector lam_new) {
......@@ -2694,15 +2697,14 @@ void Interpolate_direct::find_cell(D3vector v)
counterSamePoint++;
return;
}
badLambdaFound = false;
bool badIsBestLambda = false;
bool previousBadLambda = false;
idHexNow = -1;
vNow = v;
lambda = D3vector(-1,-1,-1);
D3vector cWSD, cESD;
D3vector cWND, cEND;
D3vector cWST, cEST;
D3vector cWNT, cENT;
D3vector boxWSD, boxENT;
......@@ -2711,58 +2713,98 @@ void Interpolate_direct::find_cell(D3vector v)
int Nx, Ny, Nz;
// prev
if (checkBox(idHexPrev,iPrev, jPrev, kPrev, v) != -1){
counterFastest++;
setPrevIndex();
//cout << " n.: " << idHexNow << " " << kNow << " " << jNow << " " << kNow << endl;
return;
if (!badLambdaFound)
{
counterFastest++;
setPrevIndex();
return;
}
else
{
previousBadLambda = true;
//std::cout << "should be in next or is at boundary!" << std::endl;
}
}
//std::cout << idHexPrev<< " " <<iPrev<< " " << jPrev<< " " << kPrev << std::endl;
if (checkBoxSurroundingOptimized(idHexPrev,iPrev, jPrev, kPrev, v) != -1)
{
if ((previousBadLambda && badLambdaFound ) || !badLambdaFound)
{
//will stay bad : probably at the boundary!
counterFast++;
setPrevIndex();
return;
}
if (checkBoxSurroundingOptimized(idHexPrev,iPrev, jPrev, kPrev, v) != -1){
counterFast++;
setPrevIndex();
//cout << " n.: " << idHexNow << " " << kNow << " " << jNow << " " << kNow << endl;
return;}
}
// prev prev
if (checkBox(idHexPrevPrev,iPrevPrev, jPrevPrev, kPrevPrev, v) != -1){
//counterFastest++;
setPrevIndex();
setPrevPrevIndex();
//cout << " n.: " << idHexNow << " " << kNow << " " << jNow << " " << kNow << endl;
counterSecondTry++;
return;}
if (checkBoxSurroundingOptimized(idHexPrevPrev,iPrevPrev, jPrevPrev, kPrevPrev, v) != -1){
//counterFast++;
setPrevIndex();
setPrevPrevIndex();
//cout << " n.: " << idHexNow << " " << kNow << " " << jNow << " " << kNow << endl;
counterSecondTry++;
return;}
if (checkBox(idHexPrevPrev,iPrevPrev, jPrevPrev, kPrevPrev, v) != -1)
{
if (!badLambdaFound)
{
setPrevIndex();
setPrevPrevIndex();
//cout << " n.: " << idHexNow << " " << kNow << " " << jNow << " " << kNow << endl;
counterSecondTry++;
return;
}
else
{
previousBadLambda = true;
//std::cout << "should be in next or is at boundary!" << std::endl;
}
}
if (checkBoxSurroundingOptimized(idHexPrevPrev,iPrevPrev, jPrevPrev, kPrevPrev, v) != -1)
{
if ((previousBadLambda && badLambdaFound ) || !badLambdaFound)
{
setPrevIndex();
setPrevPrevIndex();
//cout << " n.: " << idHexNow << " " << kNow << " " << jNow << " " << kNow << endl;
counterSecondTry++;
return;
}
}
// prev prev prev
if (checkBox(idHexPrevPrevPrev,iPrevPrevPrev, jPrevPrevPrev, kPrevPrevPrev, v) != -1){
//counterFastest++;
setPrevIndex();
setPrevPrevPrevIndex();
//cout << " n.: " << idHexNow << " " << kNow << " " << jNow << " " << kNow << endl;
counterThirdTry++;
return;}
if (checkBox(idHexPrevPrevPrev,iPrevPrevPrev, jPrevPrevPrev, kPrevPrevPrev, v) != -1)
{
if (!badLambdaFound)
{
setPrevIndex();
setPrevPrevIndex();
setPrevPrevPrevIndex();
//cout << " n.: " << idHexNow << " " << kNow << " " << jNow << " " << kNow << endl;
counterThirdTry++;
return;
}
else
{
previousBadLambda = true;
//std::cout << "should be in next or is at boundary!" << std::endl;
}
if (checkBoxSurroundingOptimized(idHexPrevPrevPrev,iPrevPrevPrev, jPrevPrevPrev, kPrevPrevPrev, v) != -1){
//counterFast++;
setPrevIndex();
setPrevPrevIndex();
setPrevPrevPrevIndex();
//cout << " n.: " << idHexNow << " " << kNow << " " << jNow << " " << kNow << endl;
counterThirdTry++;
return;}
}
if (checkBoxSurroundingOptimized(idHexPrevPrevPrev,iPrevPrevPrev, jPrevPrevPrev, kPrevPrevPrev, v) != -1)
{
if ((previousBadLambda && badLambdaFound ) || !badLambdaFound)
{
setPrevIndex();
setPrevPrevIndex();
setPrevPrevPrevIndex();
//cout << " n.: " << idHexNow << " " << kNow << " " << jNow << " " << kNow << endl;
counterThirdTry++;
return;
}
idHexNow = -1;
}
// cout << "Point looking for :";v.Print();cout<<endl;
for(int id_hex=0;id_hex<blockgrid->Give_unstructured_grid()->Give_number_hexahedra();++id_hex) {
Nx = blockgrid->Give_Nx_hexahedron(id_hex);
......@@ -2774,13 +2816,30 @@ void Interpolate_direct::find_cell(D3vector v)
for(int i=0;i<Nx;++i) {
if (0 < checkBox(id_hex,i, j, k, v))
{
if (abs(iNow - iPrev) < 2 && abs(jNow - jPrev) < 2 && abs(kNow - kPrev) < 2)
{
std::cout<<std::endl;
std::cout << "why not found?!\n";
}
int babsdbas=0;
babsdbas = 1;
// k =Nz;
// j =Ny;
// i =Nx;
// id_hex = blockgrid->Give_unstructured_grid()->Give_number_hexahedra();
}
}
}
}
if (badLambdaFound)
{
if (MIN(lambda) < 0 || MAX(lambda)>1)
{
}
else
{
assert(false && "must not happen!");
}
//lambda.Print();std::cout<<std::flush;
}
// cout << " pp: " << idHexPrevPrevPrev << " " << iPrevPrevPrev << " " << jPrevPrevPrev << " " << kPrevPrevPrev << endl;
// cout << " pp: " << idHexPrevPrev << " " << iPrevPrev << " " << jPrevPrev << " " << kPrevPrev << endl;
......@@ -2935,76 +2994,112 @@ int Interpolate_direct::checkBox(int id_Hex, int i, int j, int k, D3vector v)
cWNT = blockgrid->Give_coord_hexahedron(id_Hex,i, j+1,k+1);
cENT = blockgrid->Give_coord_hexahedron(id_Hex,i+1,j+1,k+1);
// std::cout << "delete line below " << std::endl;
// k = 1;
// i = 3;
// j = 3;
// cWSD = blockgrid->Give_coord_hexahedron(0,i, j, k );
// cESD = blockgrid->Give_coord_hexahedron(0,i+1,j ,k );
// cWND = blockgrid->Give_coord_hexahedron(0,i, j+1,k );
// cEND = blockgrid->Give_coord_hexahedron(0,i+1,j+1,k );
// cWST = blockgrid->Give_coord_hexahedron(0,i, j, k+1);
// cEST = blockgrid->Give_coord_hexahedron(0,i+1,j ,k+1);
// cWNT = blockgrid->Give_coord_hexahedron(0,i, j+1,k+1);
// cENT = blockgrid->Give_coord_hexahedron(0,i+1,j+1,k+1);
cWSD.Print();std::cout<<std::endl;
cESD.Print();std::cout<<std::endl;
cWND.Print();std::cout<<std::endl;
cEND.Print();std::cout<<std::endl;
cWST.Print();std::cout<<std::endl;
cEST.Print();std::cout<<std::endl;
cWNT.Print();std::cout<<std::endl;
cENT.Print();std::cout<<std::endl;
D3vector lam, lamTemp ;
if (trilinearInterpolationFlag)
{
D3vector lamTrilinar = trilinarInterpolation(v,id_Hex,i, j, k );
D3vector lamTrilinar = trilinarInterpolation(v,id_Hex,i, j, k );
lam = lambda_of_p_in_tet(v,cWND,cWNT,cWST,cEST);
lamTemp = lam;
if (contained_in_tet_strong(lam) && new_lam_better(lambda,lam)){typ=0;typCounter0++;}
lam = lambda_of_p_in_tet(v,cEST,cWND,cWST,cESD);
if (contained_in_tet_strong(lam) && new_lam_better(lambda,lam)){typ=1;typCounter1++;lamTemp = lam;}
if (lamTrilinar != -1)
{
typ = 6;
lambda = lamTrilinar;
typ_tet = typ;
idHexNow = id_Hex;
iNow = i;
jNow = j;
kNow = k;
return typ;
}
}
else
{
//std::cout << "\n\nstart of tet search\n\n";
lam = lambda_of_p_in_tet(v,cWND,cWNT,cWST,cEST);
//lam.Print();
//std::cout << "lambda bef :";lambda.Print();
lamTemp = lambda;
//std::cout << "new better : " << new_lam_better(lamTemp,lam) << std::endl;
//std::cout << "\nnew lam : "; lam.Print();
//std::cout << "old lam : "; lamTemp.Print();
if (contained_in_tet_strong(lam) && new_lam_better(lamTemp,lam)){typ=0;typCounter0++;lamTemp = lam;}
//std::cout << "new best : "; lamTemp.Print();
lam = lambda_of_p_in_tet(v,cEST,cWND,cWST,cESD);
//std::cout << "\nnew lam : "; lam.Print();
//std::cout << "old lam : "; lamTemp.Print();
//std::cout << "new better : " << new_lam_better(lamTemp,lam) << std::endl;
if (contained_in_tet_strong(lam) && new_lam_better(lamTemp,lam)){typ=1;typCounter1++;lamTemp = lam;}
//std::cout << "new best : "; lamTemp.Print();
lam = lambda_of_p_in_tet(v,cWND,cWSD,cWST,cESD);
//std::cout << "\nnew lam : "; lam.Print();
//std::cout << "old lam : "; lamTemp.Print();
//std::cout << "new better : " << new_lam_better(lamTemp,lam) << std::endl;
if (contained_in_tet_strong(lam) && new_lam_better(lamTemp,lam)){typ=2;typCounter2++;lamTemp = lam;}
//std::cout << "new best : "; lamTemp.Print();
lam = lambda_of_p_in_tet(v,cEST,cWND,cESD,cEND);
//std::cout << "\nnew lam : "; lam.Print();
//std::cout << "old lam : "; lamTemp.Print();
//std::cout << "new better : " << new_lam_better(lamTemp,lam) << std::endl;
if (contained_in_tet_strong(lam) && new_lam_better(lamTemp,lam)){typ=3;typCounter3++;lamTemp = lam;}
//std::cout << "new best : "; lamTemp.Print();
lam = lambda_of_p_in_tet(v,cENT,cWNT,cEST,cEND);
//std::cout << "\nnew lam : "; lam.Print();
//std::cout << "old lam : "; lamTemp.Print();
//std::cout << "new better : " << new_lam_better(lamTemp,lam) << std::endl;
if (contained_in_tet_strong(lam) && new_lam_better(lamTemp,lam)){typ=4;typCounter4++;lamTemp = lam;}
//std::cout << "new best : "; lamTemp.Print();
lam = lambda_of_p_in_tet(v,cWNT,cWND,cEST,cEND);
//std::cout << "\nnew lam : "; lam.Print();
//std::cout << "old lam : "; lamTemp.Print();
//std::cout << "new better : " << new_lam_better(lamTemp,lam) << std::endl;
if (contained_in_tet_strong(lam) && new_lam_better(lamTemp,lam)){typ=5;typCounter5++;lamTemp = lam;}
//std::cout << "new best : "; lamTemp.Print();
if( typ != -1 && new_lam_better(lambda,lamTemp))
{
if (MIN(lamTemp) < 0.0 || MAX(lamTemp) >1.0 )
{
badLambdaFound = true;
//lamTemp.Print();
}
else
{
badLambdaFound = false;
}
//lambda = lam;
// if (D3VectorNorm(lamTemp - D3vector(-0.184893695657936, 0.215428315449784, 0.533250940577171)) < 1e-5)
// {
// std::cout << "stop" << std::endl;
// }
// if (D3VectorNorm(lamTemp - D3vector(0.357354990724203, 0.152269993258782, 0.41683241494242)) < 1e-5)
// {
// std::cout << "stop" << std::endl;
// }
lambda = lamTemp;
typ_tet = typ;
idHexNow = id_Hex;
iNow = i;
jNow = j;
kNow = k;
lam = lambda_of_p_in_tet(v,cWND,cWSD,cWST,cESD);
if (contained_in_tet_strong(lam) && new_lam_better(lambda,lam)){typ=2;typCounter2++;lamTemp = lam;}
lam = lambda_of_p_in_tet(v,cEST,cWND,cESD,cEND);
if (contained_in_tet_strong(lam) && new_lam_better(lambda,lam)){typ=3;typCounter3++;lamTemp = lam;}
lam = lambda_of_p_in_tet(v,cENT,cWNT,cEST,cEND);
if (contained_in_tet_strong(lam) && new_lam_better(lambda,lam)){typ=4;typCounter4++;lamTemp = lam;}
}
}
lam = lambda_of_p_in_tet(v,cWNT,cWND,cEST,cEND);
if (contained_in_tet_strong(lam) && new_lam_better(lambda,lam)){typ=5;typCounter5++;lamTemp = lam;}
if( typ != -1 && new_lam_better(lambda,lamTemp))
{
if (MIN(lamTemp) < -0.2)
{
cout << "error !" << endl;
}
//lambda = lam;
lambda = lamTemp;
typ_tet = typ;
idHexNow = id_Hex;
iNow = i;
jNow = j;
kNow = k;
}
}
return typ;
}
......@@ -3025,28 +3120,44 @@ int Interpolate_direct::checkBoxSurrounding(int idHex, int i, int j, int k, D3ve
return -1;
}
int Interpolate_direct::checkBoxSurroundingOptimized(int idHex, int i, int j, int k, D3vector v)
int Interpolate_direct::checkBoxSurroundingOptimized(int idHex, int i, int j, int k, D3vector v, bool neglectBadLambda)
{
bool temp = badLambdaFound;
badLambdaFound = false;
int iterationPerSurrounding = 0;
//std::vector<
for(int iIter=i-1;iIter<i+2;iIter+=2){
if (checkBox(idHex,iIter, j, k, v) != -1){//cout << "iterationPerSurrounding " << iterationPerSurrounding+1<< endl;
for(int kIter=k+1;kIter>k-2;kIter-=2){
if (checkBox(idHex,i, j, kIter, v) != -1 && !badLambdaFound){//cout << "iterationPerSurrounding " << iterationPerSurrounding+1<< endl;
return 0;}
else { iterationPerSurrounding++;}
}
for(int jIter=j-1;jIter<j+2;jIter+=2){
if (checkBox(idHex,i, jIter, k, v) != -1){//cout << "iterationPerSurrounding " << iterationPerSurrounding+1<< endl;
if (badLambdaFound && !neglectBadLambda)
{
if (checkBoxSurroundingOptimized(idHexNow,iNow, jNow, kNow, v,true) != -1) return 0; else return -1;
}
for(int iIter=i-1;iIter<i+2;iIter+=2){
if (checkBox(idHex,iIter, j, k, v) != -1 && !badLambdaFound){//cout << "iterationPerSurrounding " << iterationPerSurrounding+1<< endl;
return 0;}
else { iterationPerSurrounding++;}
}
for(int kIter=k-1;kIter<k+2;kIter+=2){
if (checkBox(idHex,i, j, kIter, v) != -1){//cout << "iterationPerSurrounding " << iterationPerSurrounding+1<< endl;
if (badLambdaFound && !neglectBadLambda)
{
if (checkBoxSurroundingOptimized(idHexNow,iNow, jNow, kNow, v,true) != -1) return 0; else return -1;
}
for(int jIter=j-1;jIter<j+2;jIter+=2){
if (checkBox(idHex,i, jIter, k, v) != -1 && !badLambdaFound){//cout << "iterationPerSurrounding " << iterationPerSurrounding+1<< endl;
return 0;}
else { iterationPerSurrounding++;}
}
if (badLambdaFound && !neglectBadLambda)
{
if (checkBoxSurroundingOptimized(idHexNow,iNow, jNow, kNow, v,true) != -1) return 0; else return -1;
}
for(int kIter=k-1;kIter<k+2;kIter+=2)
for(int jIter=j-1;jIter<j+2;jIter+=2) {
if (checkBox(idHex,i, jIter, kIter, v) != -1){
if (checkBox(idHex,i, jIter, kIter, v) != -1 && !badLambdaFound){
//counterFast++;
//setPrevIndex();
......@@ -3054,9 +3165,13 @@ int Interpolate_direct::checkBoxSurroundingOptimized(int idHex, int i, int j, in
return 0;}
else { iterationPerSurrounding++;}
}
if (badLambdaFound && !neglectBadLambda)
{
if (checkBoxSurroundingOptimized(idHexNow,iNow, jNow, kNow, v,true) != -1) return 0; else return -1;
}
for(int kIter=k-1;kIter<k+2;kIter+=2)
for(int iIter=i-1;iIter<i+2;iIter+=2) {
if (checkBox(idHex,iIter, j, kIter, v) != -1){
if (checkBox(idHex,iIter, j, kIter, v) != -1 && !badLambdaFound){ //here
//counterFast++;
//setPrevIndex();
......@@ -3064,9 +3179,13 @@ int Interpolate_direct::checkBoxSurroundingOptimized(int idHex, int i, int j, in
return 0;}
else { iterationPerSurrounding++;}
}
if (badLambdaFound && !neglectBadLambda)
{
if (checkBoxSurroundingOptimized(idHexNow,iNow, jNow, kNow, v,true) != -1) return 0; else return -1;
}
for(int jIter=j-1;jIter<j+2;jIter+=2)
for(int iIter=i-1;iIter<i+2;iIter+=2) {
if (checkBox(idHex,iIter, jIter, k, v) != -1){
if (checkBox(idHex,iIter, jIter, k, v) != -1 && !badLambdaFound){
//counterFast++;
//setPrevIndex();
......@@ -3074,23 +3193,30 @@ int Interpolate_direct::checkBoxSurroundingOptimized(int idHex, int i, int j, in
return 0;}
else { iterationPerSurrounding++;}
}
if (badLambdaFound && !neglectBadLambda)
{
if (checkBoxSurroundingOptimized(idHexNow,iNow, jNow, kNow, v,true) != -1) return 0; else return -1;
}
for(int kIter=k-1;kIter<k+2;kIter+=2)
for(int jIter=j-1;jIter<j+2;jIter+=2)
for(int iIter=i-1;iIter<i+2;iIter+=2) {
if (checkBox(idHex,iIter, jIter, kIter, v) != -1){
if (checkBox(idHex,iIter, jIter, kIter, v) != -1 && !badLambdaFound){
//counterFast++;
//setPrevIndex();
//cout << "iterationPerSurrounding " << iterationPerSurrounding+1<< endl;
return 0;}
else { iterationPerSurrounding++;}
}
if (badLambdaFound && !neglectBadLambda)
{
if (checkBoxSurroundingOptimized(idHexNow,iNow, jNow, kNow, v,true) != -1) return 0; else return -1;
}
badLambdaFound = temp;
return -1;
}
......@@ -3398,23 +3524,133 @@ void Interpolate_direct::writeInterpolationEfficiency()
}
}
bool isInRange(double d0, double d1)
{
if (d0 * d1 >=0 )
{
return true;
}
return false;
}
D3vector Interpolate_direct::trilinarInterpolation(D3vector X, int id_Hex, int i, int j, int k)
{
std::cout << "to be implemented" << std::endl;
D3vector A = cEND; // 000
// cWSD = {1,1,0}; // 1,1,0 : x2
// cESD = {1,0,0}; // 1,0,0 : x1
// cWND = {0,1,0}; // 0,1,0 : x3
// cEND = {0,0,0}; // 0,0,0 : x0
// cWST = {1,1,1}; // 1,1,1 : x5
// cEST = {1,0,1}; // 1,0,1 : x6
// cWNT = {0,1,1}; // 0,1,1 : x4
// cENT = {0,0,1}; // 0,0,1 : x7
// D3vector shift{16,-26,12};
// cWSD +=shift; // 1,1,0 : x2
// cESD +=shift; // 1,0,0 : x1
// cWND +=shift; // 0,1,0 : x3
// cEND +=shift; // 0,0,0 : x0
// cWST +=shift; // 1,1,1 : x5
// cEST +=shift; // 1,0,1 : x6
// cWNT +=shift; // 0,1,1 : x4
// cENT +=shift; // 0,0,1 : x7
D3vector B = cESD - cEND; // 100 - 000
D3vector C = cWND - cEND; // 010 - 000
D3vector D = cENT - cEND; // 001 - 000
//std::cout << "to be implemented" << std::endl;
// X = (cWSD+ cESD+ cWND + cEND + cWST + cEST + cWNT + cENT ) / 8.0 ;
// X.x = X.x +15;
// X.y = X.y +15;
D3vector normalTD = cross_product(-1*B,C);
D3vector normalTDOPP = cross_product(cWST-cEST , cWNT - cWST);
//D3vector normalTDOPP2 = cross_product(cENT-cEST ,cWNT - cENT);
D3vector normalNS = cross_product(-1*C,D);
D3vector normalNSOPP = cross_product(cWST-cEST,cWST - cWSD);
//D3vector normalNSOPP2 = cross_product(cESD-cEST,cESD - cWSD);
D3vector normalEW = cross_product(-1*D,B);
D3vector normalEWOPP = cross_product(cWST-cWSD,cWST - cWNT);
//D3vector normalEWOPP2 = cross_product(cWND-cWSD,cWND - cWNT);
normalTD = normalTD / D3VectorNorm(normalTD);
normalTDOPP = normalTDOPP / D3VectorNorm(normalTDOPP);
normalNS = normalNS / D3VectorNorm(normalNS);
normalNSOPP = normalNSOPP / D3VectorNorm(normalNSOPP);
normalEW = normalEW / D3VectorNorm(normalEW);
normalEWOPP = normalEWOPP / D3VectorNorm(normalEWOPP);
double eta = 0.5;
double xi = 0.5;
double phi = 0.5;
bool fixEta = false;
bool fixXi = false;
bool fixPhi = false;
D3vector coord(eta,xi,phi);
double distTD0 = product(normalTD,cEND);
double distTD1 = product(normalTDOPP,cWST);
double distTDMid = distTD0-product(normalTD,X);
double distTDMid2 = distTD1-product(normalTDOPP,X);
double distNS0 = product(normalNS,cEND);
double distNS1 = product(normalNSOPP,cWST);
double distNSMid = distNS0-product(normalNS,X);
double distNSMid2 = distNS1 -product(normalNSOPP,X);
double distEW0 = product(normalEW,cEND);
double distEW1 = product(normalEWOPP,cWST);
double distEWMid = distEW0-product(normalEW,X);
double distEWMid2 = distEW1 -product(normalEWOPP,X);
// std::cout << "in range " << isInRange(distTDMid,distTDMid2) << std::endl;