Commit 4f3e9058 authored by Phillip Lino Rall's avatar Phillip Lino Rall
Browse files

trilinear interpolation for grids available, but doesn't always give good...

trilinear interpolation for grids available, but doesn't always give good results. needs to be improved
parent 1178403f
...@@ -372,6 +372,11 @@ inline std::vector<int> find_p_quad ( dir3D quad, Hexahedron_el *hex ) ...@@ -372,6 +372,11 @@ inline std::vector<int> find_p_quad ( dir3D quad, Hexahedron_el *hex )
// WDed, EDed, WTed, ETed, // WDed, EDed, WTed, ETed,
// SWed, SEed, NWed, NEed }; // SWed, SEed, NWed, NEed };
std::vector<int> ijk = {1 , 1 , 1}; std::vector<int> ijk = {1 , 1 , 1};
// for (int iter = 0 ; iter < 12; iter++)
// {
// std::cout << hex->orientation[iter] << "\n" ;
// }
// std::cout << std::endl;
if (quad == Wdir3D) if (quad == Wdir3D)
{ {
...@@ -550,14 +555,45 @@ D3vector Blockgrid::Give_coord_hexahedron ( int id_hex, ...@@ -550,14 +555,45 @@ D3vector Blockgrid::Give_coord_hexahedron ( int id_hex,
int quadWaddK = (p_quad_Wdir3D[2] == 1) ? 0 : nEdgeZ; int quadWaddK = (p_quad_Wdir3D[2] == 1) ? 0 : nEdgeZ;
// std::cout <<"Nx " << Give_Nx_quadrangle ( idW ) << " Ny " << Give_Ny_quadrangle ( idW ) << "\n";
// std::cout <<"Nx " << Give_Nx_quadrangle ( idE ) << " Ny " << Give_Ny_quadrangle ( idE ) << "\n";
// std::cout <<"Nx " << Give_Nx_quadrangle ( idS ) << " Ny " << Give_Ny_quadrangle ( idS ) << "\n";
// std::cout <<"Nx " << Give_Nx_quadrangle ( idN ) << " Ny " << Give_Ny_quadrangle ( idN ) << "\n";
// std::cout <<"Nx " << Give_Nx_quadrangle ( idD ) << " Ny " << Give_Ny_quadrangle ( idD ) << "\n";
// std::cout <<"Nx " << Give_Nx_quadrangle ( idT ) << " Ny " << Give_Ny_quadrangle ( idT ) << std::endl;
// std::cout << std::endl;
// std::cout << " transface coordiantes: "<< std::endl;
//with modification //with modification
QPsi[0] = Give_coord_quadrangle( idW , quadWaddJ + j*p_quad_Wdir3D[1] , quadWaddK + k*p_quad_Wdir3D[2]);
QPsi[1] = Give_coord_quadrangle( idE , quadEaddJ + j*p_quad_Edir3D[1] , quadEaddK + k*p_quad_Edir3D[2]);
QPsi[2] = Give_coord_quadrangle( idS , quadSaddI + i*p_quad_Sdir3D[0] , quadSaddK + k*p_quad_Sdir3D[2]);
QPsi[3] = Give_coord_quadrangle( idN , quadNaddI + i*p_quad_Ndir3D[0] , quadNaddK + k*p_quad_Ndir3D[2]);
QPsi[4] = Give_coord_quadrangle( idD , quadDaddI + i*p_quad_Ddir3D[0] , quadDaddJ + j*p_quad_Ddir3D[1]);
QPsi[5] = Give_coord_quadrangle( idT , quadTaddI + i*p_quad_Tdir3D[0] , quadTaddJ + j*p_quad_Tdir3D[1]);
bool invertTD = false;
bool invertNS = false;
bool invertEW = false;
if (!(nEdgeY == Give_Nx_quadrangle ( idW ) && nEdgeZ == Give_Ny_quadrangle ( idW )))
{
invertEW = true;
}
if (!(nEdgeX == Give_Nx_quadrangle ( idS ) && nEdgeZ == Give_Ny_quadrangle ( idS ) ))
{
invertNS = true;
}
if (!(nEdgeX == Give_Nx_quadrangle ( idD ) && nEdgeY == Give_Ny_quadrangle ( idD )))
{
invertTD = true;
}
QPsi[0] = Give_coord_quadrangle( idW , quadWaddJ + j*p_quad_Wdir3D[1] , quadWaddK + k*p_quad_Wdir3D[2],invertEW);
QPsi[1] = Give_coord_quadrangle( idE , quadEaddJ + j*p_quad_Edir3D[1] , quadEaddK + k*p_quad_Edir3D[2],invertEW);
QPsi[2] = Give_coord_quadrangle( idS , quadSaddI + i*p_quad_Sdir3D[0] , quadSaddK + k*p_quad_Sdir3D[2],invertNS);
QPsi[3] = Give_coord_quadrangle( idN , quadNaddI + i*p_quad_Ndir3D[0] , quadNaddK + k*p_quad_Ndir3D[2],invertNS);
QPsi[4] = Give_coord_quadrangle( idD , quadDaddI + i*p_quad_Ddir3D[0] , quadDaddJ + j*p_quad_Ddir3D[1],invertTD);
QPsi[5] = Give_coord_quadrangle( idT , quadTaddI + i*p_quad_Tdir3D[0] , quadTaddJ + j*p_quad_Tdir3D[1],invertTD);
//delete again!!!
// QPsi[5] = Give_coord_quadrangle( idT , quadTaddJ + j*p_quad_Tdir3D[1] , quadTaddI + i*p_quad_Tdir3D[0]);
//std::cout << std::endl;
// Version: 1 x Flaechen + 2 x Kanten + 3 x Ecken // Version: 1 x Flaechen + 2 x Kanten + 3 x Ecken
// ///////////////////////////////////////////////////////// // /////////////////////////////////////////////////////////
Pres = (QPsi[0] + eta * (QPsi[1]-QPsi[0]) + QPsi[2] + xi * (QPsi[3]-QPsi[2]) + QPsi[4] + phi * (QPsi[5]-QPsi[4])); Pres = (QPsi[0] + eta * (QPsi[1]-QPsi[0]) + QPsi[2] + xi * (QPsi[3]-QPsi[2]) + QPsi[4] + phi * (QPsi[5]-QPsi[4]));
...@@ -597,21 +633,21 @@ D3vector Blockgrid::Give_coord_hexahedron ( int id_hex, ...@@ -597,21 +633,21 @@ D3vector Blockgrid::Give_coord_hexahedron ( int id_hex,
if ( md==0 ) // EW if ( md==0 ) // EW
{ {
PSW = Give_coord_quadrangle(idS,quadSaddI+i*p_quad_Sdir3D[0] ,quadSaddK); PSW = Give_coord_quadrangle(idS,quadSaddI+i*p_quad_Sdir3D[0] ,quadSaddK, invertNS);
if(coordianteFromCorners) if(coordianteFromCorners)
PSW = Psi[SDed]; PSW = Psi[SDed];
PSE = Give_coord_quadrangle(idN,quadNaddI+i * p_quad_Ndir3D[0],quadNaddK); PSE = Give_coord_quadrangle(idN,quadNaddI+i * p_quad_Ndir3D[0],quadNaddK, invertNS);
if(coordianteFromCorners) if(coordianteFromCorners)
PSE = Psi[NDed]; PSE = Psi[NDed];
PNW = Give_coord_quadrangle(idS,quadSaddI+i*p_quad_Sdir3D[0] ,nEdgeZ-quadSaddK); PNW = Give_coord_quadrangle(idS,quadSaddI+i*p_quad_Sdir3D[0] ,nEdgeZ-quadSaddK, invertNS);
if(coordianteFromCorners) if(coordianteFromCorners)
PNW = Psi[STed]; PNW = Psi[STed];
PNE = Give_coord_quadrangle(idN,quadNaddI+i * p_quad_Ndir3D[0],nEdgeZ-quadNaddK); PNE = Give_coord_quadrangle(idN,quadNaddI+i * p_quad_Ndir3D[0],nEdgeZ-quadNaddK, invertNS);
if(coordianteFromCorners) if(coordianteFromCorners)
PNE = Psi[NTed]; PNE = Psi[NTed];
...@@ -625,21 +661,21 @@ D3vector Blockgrid::Give_coord_hexahedron ( int id_hex, ...@@ -625,21 +661,21 @@ D3vector Blockgrid::Give_coord_hexahedron ( int id_hex,
{ {
PSW = Give_coord_quadrangle(idD,quadDaddI,quadDaddJ+j*p_quad_Ddir3D[1]); PSW = Give_coord_quadrangle(idD,quadDaddI,quadDaddJ+j*p_quad_Ddir3D[1], invertTD);
if(coordianteFromCorners) if(coordianteFromCorners)
PSW = Psi[WDed]; PSW = Psi[WDed];
PSE = Give_coord_quadrangle(idD,nEdgeX-quadDaddI,quadDaddJ+j*p_quad_Ddir3D[1]); PSE = Give_coord_quadrangle(idD,nEdgeX-quadDaddI,quadDaddJ+j*p_quad_Ddir3D[1], invertTD);
if(coordianteFromCorners) if(coordianteFromCorners)
PSE = Psi[EDed]; PSE = Psi[EDed];
PNW = Give_coord_quadrangle(idT,quadTaddI,quadTaddJ+j*p_quad_Tdir3D[1]); PNW = Give_coord_quadrangle(idT,quadTaddI,quadTaddJ+j*p_quad_Tdir3D[1], invertTD);
if(coordianteFromCorners) if(coordianteFromCorners)
PNW = Psi[WTed]; PNW = Psi[WTed];
PNE = Give_coord_quadrangle(idT,nEdgeX-quadTaddI,quadTaddJ+j*p_quad_Tdir3D[1]); PNE = Give_coord_quadrangle(idT,nEdgeX-quadTaddI,quadTaddJ+j*p_quad_Tdir3D[1], invertTD);
if(coordianteFromCorners) if(coordianteFromCorners)
PNE = Psi[ETed]; PNE = Psi[ETed];
...@@ -654,20 +690,20 @@ D3vector Blockgrid::Give_coord_hexahedron ( int id_hex, ...@@ -654,20 +690,20 @@ D3vector Blockgrid::Give_coord_hexahedron ( int id_hex,
{ {
PSW = Give_coord_quadrangle(idW,quadWaddJ,quadWaddK+k*p_quad_Wdir3D[2]); PSW = Give_coord_quadrangle(idW,quadWaddJ,quadWaddK+k*p_quad_Wdir3D[2], invertEW);
if(coordianteFromCorners) if(coordianteFromCorners)
PSW = Psi[SWed]; PSW = Psi[SWed];
PSE = Give_coord_quadrangle(idE,quadEaddJ,quadEaddK+k*p_quad_Edir3D[2]); PSE = Give_coord_quadrangle(idE,quadEaddJ,quadEaddK+k*p_quad_Edir3D[2], invertEW);
if(coordianteFromCorners) if(coordianteFromCorners)
PSE = Psi[SEed]; PSE = Psi[SEed];
PNW = Give_coord_quadrangle(idW,nEdgeY-quadWaddJ,quadWaddK+k*p_quad_Wdir3D[2]); PNW = Give_coord_quadrangle(idW,nEdgeY-quadWaddJ,quadWaddK+k*p_quad_Wdir3D[2], invertEW);
if(coordianteFromCorners) if(coordianteFromCorners)
PNW = Psi[NWed]; PNW = Psi[NWed];
PNE = Give_coord_quadrangle(idE,nEdgeY-quadEaddJ,quadEaddK+k*p_quad_Edir3D[2]); PNE = Give_coord_quadrangle(idE,nEdgeY-quadEaddJ,quadEaddK+k*p_quad_Edir3D[2], invertEW);
if(coordianteFromCorners) if(coordianteFromCorners)
PNE = Psi[NEed]; PNE = Psi[NEed];
...@@ -764,7 +800,7 @@ D3vector Blockgrid::Give_coord_hexahedron ( int id_hex, ...@@ -764,7 +800,7 @@ D3vector Blockgrid::Give_coord_hexahedron ( int id_hex,
return Pres + hex->transform_hex ( id_hex,eta,xi,phi ); return Pres + hex->transform_hex ( id_hex,eta,xi,phi );
} }
D3vector Blockgrid::Give_coord_quadrangle ( int id, int i, int j ) const D3vector Blockgrid::Give_coord_quadrangle (int id, int i, int j, bool invert) const
{ {
if(bg_coord != NULL) if(bg_coord->blockgrid_quad_coordinates_calculated) if(bg_coord != NULL) if(bg_coord->blockgrid_quad_coordinates_calculated)
{ {
...@@ -775,6 +811,8 @@ D3vector Blockgrid::Give_coord_quadrangle ( int id, int i, int j ) const ...@@ -775,6 +811,8 @@ D3vector Blockgrid::Give_coord_quadrangle ( int id, int i, int j ) const
} }
double eta, xi; double eta, xi;
D3vector PSW, PSE, PNW, PNE, PW, PE, PN, PS; D3vector PSW, PSE, PNW, PNE, PW, PE, PN, PS;
D3vector PTransFace, PWithoutTrans; D3vector PTransFace, PWithoutTrans;
...@@ -783,9 +821,30 @@ D3vector Blockgrid::Give_coord_quadrangle ( int id, int i, int j ) const ...@@ -783,9 +821,30 @@ D3vector Blockgrid::Give_coord_quadrangle ( int id, int i, int j ) const
quad = this->ug->Give_quadrangle ( id ); quad = this->ug->Give_quadrangle ( id );
double Nx = (double) Give_Nx_quadrangle ( id );
double Ny = ( double ) Give_Ny_quadrangle ( id );
// std::cout << "nx ny " << Nx << " " << Ny << std::endl;
// if (id == 23 && Nx != Ny)
// {
// std::cout << "break";
// }
// std::cout << "nx ny " << Nx << " " << Ny << std::endl;
if (invert)
{
double temp = Nx; Nx = Ny; Ny = temp;
}
//std::cout << "nx ny " << Nx << " " << Ny << std::endl;
eta = ( double ) i/ Nx;
xi = ( double ) j/ Ny;
if (eta > 1 || xi > 1)
{
std::cout << "debug : flag is " << invert << ", but should be " << !invert << std::endl;
std::cout << std::endl;
}
eta = ( double ) i/ ( double ) Give_Nx_quadrangle ( id );
xi = ( double ) j/ ( double ) Give_Ny_quadrangle ( id );
// @todo : so sollte es sein // @todo : so sollte es sein
PSW = quad->Give_coord ( SWdir2D ); PSW = quad->Give_coord ( SWdir2D );
...@@ -814,7 +873,9 @@ D3vector Blockgrid::Give_coord_quadrangle ( int id, int i, int j ) const ...@@ -814,7 +873,9 @@ D3vector Blockgrid::Give_coord_quadrangle ( int id, int i, int j ) const
if (quad->transform != NULL) if (quad->transform != NULL)
{ {
PTransFace = quad->transform ( eta, xi, quad->shiftPointer(ug->Get_pointer_global_data())); PTransFace = quad->transform ( eta, xi, quad->shiftPointer(ug->Get_pointer_global_data()));
PTransFace = PSW + (PSE-PSW) * eta + (PNW-PSW)*xi + (PNE-PSE-PNW+PSW)*xi*eta + PTransFace; //std::cout << "PTransFace ";PTransFace.Print();
D3vector TEMP = PSW + (PSE-PSW) * eta + (PNW-PSW)*xi + (PNE-PSE-PNW+PSW)*xi*eta ;
PTransFace = TEMP + PTransFace;// + PTransFace;
return PTransFace; return PTransFace;
} }
...@@ -849,7 +910,6 @@ D3vector Blockgrid::Give_coord_edge ( int id, int i ) const ...@@ -849,7 +910,6 @@ D3vector Blockgrid::Give_coord_edge ( int id, int i ) const
if(bg_coord != NULL) if(bg_coord->blockgrid_edge_coordinates_calculated) if(bg_coord != NULL) if(bg_coord->blockgrid_edge_coordinates_calculated)
{ {
return bg_coord->blockgrid_edge_coordinates.at(id).at(i) ; return bg_coord->blockgrid_edge_coordinates.at(id).at(i) ;
} }
...@@ -903,49 +963,85 @@ D3vector Blockgrid::Give_coord_edge ( int id, int i ) const ...@@ -903,49 +963,85 @@ D3vector Blockgrid::Give_coord_edge ( int id, int i ) const
} }
} }
} }
ug->Set_edge_to_hex_id(id, idHex);
ug->Set_edge_to_quad_id(id,quadId); ug->Set_edge_to_quad_id(id,quadId);
ug->Set_edge_to_quad_dir(id,(int)dirForQuad); ug->Set_edge_to_quad_dir(id,(int)dirForQuad);
} }
else else
{ {
idHex = ug->Give_edge_to_hex_id(id);
quadId = ug->Give_edge_to_quad_id(id); quadId = ug->Give_edge_to_quad_id(id);
dirForQuad = (dir3D)ug->Give_edge_to_quad_dir(id); dirForQuad = (dir3D)ug->Give_edge_to_quad_dir(id);
quad = ug->Give_quadrangle(quadId); quad = ug->Give_quadrangle(quadId);
} }
int Nx = Give_Nx_quadrangle(dirForQuad); int Nx = Give_Nx_quadrangle(quadId);
int Ny = Give_Ny_quadrangle(dirForQuad); int Ny = Give_Ny_quadrangle(quadId);
int nEdgeX = Give_Nx_hexahedron ( idHex );
int nEdgeY = Give_Ny_hexahedron ( idHex );
int nEdgeZ = Give_Nz_hexahedron ( idHex );
bool invert = false;
// std::cout << "dirForQuad " << dirForQuad << std::endl;
// std::cout << "Nx " << Nx << " Ny " << Ny << std::endl;
// if (!(Nx == Give_Nx_quadrangle ( dirForQuad ) && Ny == Give_Ny_quadrangle ( dirForQuad ) ))
// {
// invert = true;
// }
// if (dirForQuad == Edir3D || dirForQuad == Wdir3D)
// {
// if (!(nEdgeY == Give_Nx_quadrangle ( dirForQuad ) && nEdgeZ == Give_Ny_quadrangle ( dirForQuad ) ))
// {
// invert = true;
// }
// }
// else if (dirForQuad == Ndir3D || dirForQuad == Sdir3D)
// {
// if (!(nEdgeX == Give_Nx_quadrangle ( dirForQuad ) && nEdgeZ == Give_Ny_quadrangle ( dirForQuad ) ))
// {
// invert = true;
// }
// }
// else if (dirForQuad == Tdir3D || dirForQuad == Ddir3D)
// {
// if (!(nEdgeX == Give_Nx_quadrangle ( dirForQuad ) && nEdgeY == Give_Ny_quadrangle ( dirForQuad )))
// {
// invert = true;
// }
// }
// else
// {
// std::cout << "should not happen, dirForQuad only has 6 possibilities\n";
// }
if (Nx == nEdgeX && Ny == nEdgeY || Nx == nEdgeX && Ny == nEdgeZ || Nx == nEdgeY && Ny == nEdgeZ )
{
invert = false;
}
else
{
invert = true;
}
invert = false;
if (id == quad->Give_id_edge(Wdir2D)) //Wdir2D if (id == quad->Give_id_edge(Wdir2D)) //Wdir2D
{ {
Psi = Give_coord_quadrangle(quadId,0,i,invert);
Psi = Give_coord_quadrangle(quadId,0,i);
} }
if (id == quad->Give_id_edge(Edir2D)) //Edir2D if (id == quad->Give_id_edge(Edir2D)) //Edir2D
{ {
Psi = Give_coord_quadrangle(quadId,Nx,i,invert);
Psi = Give_coord_quadrangle(quadId,Nx,i);
} }
if (id == quad->Give_id_edge(Sdir2D)) //Sdir2D if (id == quad->Give_id_edge(Sdir2D)) //Sdir2D
{ {
Psi = Give_coord_quadrangle(quadId,i,0,invert);
Psi = Give_coord_quadrangle(quadId,i,0);
} }
if (id == quad->Give_id_edge(Ndir2D)) //Ndir2D if (id == quad->Give_id_edge(Ndir2D)) //Ndir2D
{ {
Psi = Give_coord_quadrangle(quadId,i,Ny,invert);
Psi = Give_coord_quadrangle(quadId,i,Ny);
} }
return Psi; return Psi;
...@@ -1066,6 +1162,8 @@ void Blockgrid_coordinates::init_blockgrid_coordinates() ...@@ -1066,6 +1162,8 @@ void Blockgrid_coordinates::init_blockgrid_coordinates()
} }
} }
//test:: ??? Nx = 3????
blockgrid_edge_coordinates.resize(bg->Give_unstructured_grid()->Give_number_edges()); blockgrid_edge_coordinates.resize(bg->Give_unstructured_grid()->Give_number_edges());
for (int id_edge = 0 ; id_edge < bg->Give_unstructured_grid()->Give_number_edges() ; id_edge++) for (int id_edge = 0 ; id_edge < bg->Give_unstructured_grid()->Give_number_edges() ; id_edge++)
{ {
......
...@@ -114,7 +114,7 @@ class Blockgrid { ...@@ -114,7 +114,7 @@ class Blockgrid {
void Set_blockgrid_coordinates(Blockgrid_coordinates* bg_coord_ ) { bg_coord = bg_coord_; }; void Set_blockgrid_coordinates(Blockgrid_coordinates* bg_coord_ ) { bg_coord = bg_coord_; };
D3vector Give_coord_hexahedron(int id, int i, int j, int k) const ; D3vector Give_coord_hexahedron(int id, int i, int j, int k) const ;
D3vector Give_coord_quadrangle(int id, int i, int j) const; D3vector Give_coord_quadrangle(int id, int i, int j, bool invert = false) const;
D3vector Give_coord_edge(int id, int i) const; D3vector Give_coord_edge(int id, int i) const;
D3vector Give_coord_point(int id) const; D3vector Give_coord_point(int id) const;
......
...@@ -450,6 +450,7 @@ D3vector transform_left_lens_diag_NE_quad ( double t1, double t2, double* pointe ...@@ -450,6 +450,7 @@ D3vector transform_left_lens_diag_NE_quad ( double t1, double t2, double* pointe
} }
D3vector transform_diag_inner_faces_NE_quad( double t1, double t2, double* pointer_global_data) { D3vector transform_diag_inner_faces_NE_quad( double t1, double t2, double* pointer_global_data) {
return D3vector{0,0,0};
double R_global_data = pointer_global_data[0]; double R_global_data = pointer_global_data[0];
double r_global_data = pointer_global_data[1]; double r_global_data = pointer_global_data[1];
double curvatureLeft_global_data = pointer_global_data[2]; double curvatureLeft_global_data = pointer_global_data[2];
...@@ -507,6 +508,7 @@ D3vector transform_diag_inner_faces_NE_quad( double t1, double t2, double* point ...@@ -507,6 +508,7 @@ D3vector transform_diag_inner_faces_NE_quad( double t1, double t2, double* point
} }
D3vector transform_diag_inner_faces_NE_quad_cut( double t1, double t2, double* pointer_global_data) { D3vector transform_diag_inner_faces_NE_quad_cut( double t1, double t2, double* pointer_global_data) {
return D3vector{0,0,0};
double R_global_data = pointer_global_data[0]; double R_global_data = pointer_global_data[0];
double r_global_data = pointer_global_data[1]; double r_global_data = pointer_global_data[1];
double curvatureLeft_global_data = pointer_global_data[2]; double curvatureLeft_global_data = pointer_global_data[2];
...@@ -562,6 +564,7 @@ D3vector transform_diag_inner_faces_NE_quad_cut( double t1, double t2, double* p ...@@ -562,6 +564,7 @@ D3vector transform_diag_inner_faces_NE_quad_cut( double t1, double t2, double* p
zLeft = 0; zLeft = 0;
} }
double z = zRight * (t1) + zLeft * (1-t1) ; double z = zRight * (t1) + zLeft * (1-t1) ;
if (( t1 == 0 && t2 == 0)||( t1 == 1 && t2 == 0)||( t1 == 0 && t2 == 1)||( t1 == 1 && t2 == 1)) if (( t1 == 0 && t2 == 0)||( t1 == 1 && t2 == 0)||( t1 == 0 && t2 == 1)||( t1 == 1 && t2 == 1))
{ {
...@@ -658,8 +661,12 @@ D3vector transform_right_lens_diag_NE_quad ( double t1, double t2, double* point ...@@ -658,8 +661,12 @@ D3vector transform_right_lens_diag_NE_quad ( double t1, double t2, double* point
radiusSquared) radiusSquared)
+ sign*(-curvatureRight_global_data +( (1-t)*(thickness_global_data-z_right_inner_global_data) + t * (thickness_global_data-z_right_outer_global_data))))); + sign*(-curvatureRight_global_data +( (1-t)*(thickness_global_data-z_right_inner_global_data) + t * (thickness_global_data-z_right_outer_global_data)))));
//std::cout << "z inner " << z << std::endl;
if (std::isnan(x) || std::isnan(y) || std::isnan(z)) if (std::isnan(x) || std::isnan(y) || std::isnan(z))
{ z = 0;} { z = 0;}
// std::cout << "transform_right_lens_diag_NE_quad ";
// D3vector ( x,y,z).Print();
// std::cout << std::endl;
return D3vector ( x,y,z); return D3vector ( x,y,z);
} }
...@@ -1060,13 +1067,15 @@ D3vector transform_left_lens_diag_NE_quad_cut ( double t1, double t2, double* po ...@@ -1060,13 +1067,15 @@ D3vector transform_left_lens_diag_NE_quad_cut ( double t1, double t2, double* po
x = x*t2; x = x*t2;
double y = R_global_data * (cos ( t1*0.5*M_PI )- ( 1-t1 )) ; double y = R_global_data * (cos ( t1*0.5*M_PI )- ( 1-t1 )) ;
y = y*t2; y = y*t2;
double xT = x;
double yT = y;
//added due to bended inner part : //added due to bended inner part :
double xAdd = r_global_data * ( sin ( t1*0.5*M_PI )-t1); double xAdd = r_global_data * ( sin ( t1*0.5*M_PI )-t1);
double yAdd = r_global_data *(cos ( t1*0.5*M_PI )- ( 1-t1 )) ; double yAdd = r_global_data *(cos ( t1*0.5*M_PI )- ( 1-t1 )) ;
x += (1-t2) * xAdd; x += (1-t2) * xAdd;
y += (1-t2) * yAdd; y += (1-t2) * yAdd;
double actualX = 0 + r_global_data * t1 + 0 * t2 + (R_global_data-r_global_data) * t1 * t2; double actualX = 0 + r_global_data * t1 + 0 * t2 + (R_global_data-r_global_data) * t1 * t2;
double actualY = r_global_data - r_global_data * t1 + (R_global_data-r_global_data) * t2 - (R_global_data-r_global_data) * t1 * t2; double actualY = r_global_data - r_global_data * t1 + (R_global_data-r_global_data) * t2 - (R_global_data-r_global_data) * t1 * t2;
...@@ -1079,12 +1088,14 @@ D3vector transform_left_lens_diag_NE_quad_cut ( double t1, double t2, double* po ...@@ -1079,12 +1088,14 @@ D3vector transform_left_lens_diag_NE_quad_cut ( double t1, double t2, double* po
if (z_left_inner_global_data == z_left_outer_global_data) if (z_left_inner_global_data == z_left_outer_global_data)
{ {
//std::cout << "set to zero again!\n";
z = 0; z = 0;
} }
if (std::isnan(z)) if (std::isnan(z))
{z = 0;} {z = 0;}
//return D3vector ( xT,yT,z);
return D3vector ( x,y,z); return D3vector ( x,y,z);
} }
...@@ -1121,6 +1132,8 @@ D3vector transform_right_lens_diag_NE_quad_cut ( double t1, double t2, double* p ...@@ -1121,6 +1132,8 @@ D3vector transform_right_lens_diag_NE_quad_cut ( double t1, double t2, double* p
double y = R_global_data * (cos ( t1*0.5*M_PI )- ( 1-t1 )) ; double y = R_global_data * (cos ( t1*0.5*M_PI )- ( 1-t1 )) ;
y = y*t2; y = y*t2;
double xT = x;
double yT = y;
//added due to bended inner part : //added due to bended inner part :
double xAdd = r_global_data * ( sin ( t1*0.5*M_PI )-t1); double xAdd = r_global_data * ( sin ( t1*0.5*M_PI )-t1);
double yAdd = r_global_data *(cos ( t1*0.5*M_PI )- ( 1-t1 )) ; double yAdd = r_global_data *(cos ( t1*0.5*M_PI )- ( 1-t1 )) ;
...@@ -1133,11 +1146,18 @@ D3vector transform_right_lens_diag_NE_quad_cut ( double t1, double t2, double* p ...@@ -1133,11 +1146,18 @@ D3vector transform_right_lens_diag_NE_quad_cut ( double t1, double t2, double* p
double radiusSquared = ((x+actualX)*(x+actualX)+(y+actualY)*(y+actualY)); double radiusSquared = ((x+actualX)*(x+actualX)+(y+actualY)*(y+actualY));
// std::cout << "(-curvatureRight_global_data +( (1-t)*(thickness_global_data-z_right_inner_global_data) + t * (thickness_global_data-z_right_outer_global_data))) " << (-curvatureRight_global_data +( (1-t)*(thickness_global_data-z_right_inner_global_data) + t * (thickness_global_data-z_right_outer_global_data))) << "\n";
// std::cout << "t1,2 :: " << t1 << " " << t2 << "\n";
//std::cout << "radius " << sqrt(radiusSquared) << std::endl;
double z = offsetZ_global_data+ double z = offsetZ_global_data+
sign*( ( sqrt(pow(curvatureRight_global_data,2)- sign*( ( sqrt(pow(curvatureRight_global_data,2)-
radiusSquared) radiusSquared)
+ sign*(-curvatureRight_global_data +( (1-t)*(thickness_global_data-z_right_inner_global_data) + t * (thickness_global_data-z_right_outer_global_data))))); + sign*(-curvatureRight_global_data +( (1-t)*(thickness_global_data-z_right_inner_global_data) + t * (thickness_global_data-z_right_outer_global_data)))));
// z =radiusSquared;
//z *= 10;
//z = 1;
//std::cout << "z outer" << z << std::endl;
// double ADAD = offsetZ_global_data+sign*(-1*( sqrt(pow(curvatureLeft_global_data,2)-pow(r_global_data * (1-t) + R_global_data * t,2)) +sign *(-curvatureLeft_global_data +( (1-t)*z_left_inner_global_data + t * z_left_outer_global_data))));
if (z_right_inner_global_data == z_right_outer_global_data) if (z_right_inner_global_data == z_right_outer_global_data)
{ {
...@@ -1146,7 +1166,19 @@ D3vector transform_right_lens_diag_NE_quad_cut ( double t1, double t2, double* p ...@@ -1146,7 +1166,19 @@ D3vector transform_right_lens_diag_NE_quad_cut ( double t1, double t2, double* p
if (std::isnan(x) || std::isnan(y) || std::isnan(z)) if (std::isnan(x) || std::isnan(y) || std::isnan(z))
{ z = 0;} { z = 0;}
// std::cout << "transform_