#include "../mympi.h"
#include "../abbrevi.h"
#include "../parameter.h"
#include "../math_lib/math_lib.h"
#include "../basics/basic.h"
#include "elements.h"
#include "parti.h"
#include "ug.h"
#include <algorithm>
#include "examples_ug_optics.h"

////////////////////////////////////////////////////////////////////
// cylinder : params already defined, but for understanding the code (me, phillip) this stays here.
////////////////////////////////////////////////////////////////////

D3vector transform_quadrangle_NULL( double t1,double t2, double* pointer_global_data ) {
  return D3vector (0,0,0);
}

D3vector transform_outer_boundary_SE( double t1,double t2, double* pointer_global_data ) {

    double R_global_data = pointer_global_data[0];
    double x = (cos ( t1*0.5*M_PI )- ( 1-t1 )) ;
    double y = -( sin ( t1*0.5*M_PI )-t1);
    double z = 0;
    return R_global_data*D3vector ( x,y,z);
}
D3vector transform_outer_boundary_SW( double t1,double t2, double* pointer_global_data ) {

    double R_global_data = pointer_global_data[0];
    double x = -(cos ( t1*0.5*M_PI )- ( 1-t1 )) ;
    double y = -( sin ( t1*0.5*M_PI )-t1);
    double z = 0;
    return R_global_data*D3vector ( x,y,z);
}

D3vector transform_outer_boundary_NE( double t1,double t2, double* pointer_global_data ) {

    double R_global_data = pointer_global_data[0];
    double x = (cos ( t1*0.5*M_PI )- ( 1-t1 )) ;
    double y = ( sin ( t1*0.5*M_PI )-t1);
    double z = 0;
    return R_global_data*D3vector ( x,y,z);
}
D3vector transform_outer_boundary_NW( double t1,double t2, double* pointer_global_data ) {

    double R_global_data = pointer_global_data[0];
    double x = -( sin ( t1*0.5*M_PI )-t1);
    double y = (cos ( t1*0.5*M_PI )- ( 1-t1 )) ;
    double z = 0;
    return R_global_data*D3vector ( x,y,z);
}

D3vector transform_outer_boundary_SE_cut( double t1,double t2, double* pointer_global_data ) {

    double R_outer = pointer_global_data[12];
    double x = (cos ( t1*0.5*M_PI )- ( 1-t1 )) ;
    double y = -( sin ( t1*0.5*M_PI )-t1);
    double z = 0;
    return R_outer*D3vector ( x,y,z);
}
D3vector transform_outer_boundary_SW_cut( double t1,double t2, double* pointer_global_data ) {

    double R_outer = pointer_global_data[12];
    double x = -(cos ( t1*0.5*M_PI )- ( 1-t1 )) ;
    double y = -( sin ( t1*0.5*M_PI )-t1);
    double z = 0;
    return R_outer*D3vector ( x,y,z);
}

D3vector transform_outer_boundary_NE_cut( double t1,double t2, double* pointer_global_data ) {

    double R_outer = pointer_global_data[12];
    double x = (cos ( t1*0.5*M_PI )- ( 1-t1 )) ;
    double y = ( sin ( t1*0.5*M_PI )-t1);
    double z = 0;
    return R_outer*D3vector ( x,y,z);
}
D3vector transform_outer_boundary_NW_cut( double t1,double t2, double* pointer_global_data ) {

    double R_outer = pointer_global_data[12];
    double x = -( sin ( t1*0.5*M_PI )-t1);
    double y = (cos ( t1*0.5*M_PI )- ( 1-t1 )) ;
    double z = 0;
    return R_outer*D3vector ( x,y,z);
}

D3vector transform_right_lens_diag_NW_quad ( double t1, double t2, double* pointer_global_data) {

    double R_global_data              = pointer_global_data[0];
    double r_global_data              = pointer_global_data[1];
    double curvatureLeft_global_data  = pointer_global_data[2];
    double curvatureRight_global_data = pointer_global_data[3];
    double thickness_global_data      = pointer_global_data[4];
    double offsetZ_global_data        = pointer_global_data[5];
    double z_left_inner_global_data   = pointer_global_data[6];
    double z_left_outer_global_data   = pointer_global_data[7];
    double z_right_inner_global_data  = pointer_global_data[8];
    double z_right_outer_global_data  = pointer_global_data[9];
    double cut_edge_from_left         = pointer_global_data[10];
    double cut_edge_from_right        = pointer_global_data[11];

    double temp = t1;
    t1 = t2;
    t2 = temp;
    double t = t2;
   // t2 = (1-t2);
    double x = -R_global_data*( sin ( t1*0.5*M_PI )-t1);
    x = x*t2;
    double y = R_global_data*(cos ( t1*0.5*M_PI )- ( 1-t1 )) ;
    y = y*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;

    actualX = -actualX;
    double radiusSquared = ((x+actualX)*(x+actualX)+(y+actualY)*(y+actualY));
    double sign = (curvatureRight_global_data > 0) ? 1 : ((curvatureRight_global_data < 0) ? -1 : 0)  ;

    double z = offsetZ_global_data+
                sign*(  ( sqrt(pow(curvatureRight_global_data,2)-
                               radiusSquared)
                          + sign*(-curvatureRight_global_data +( (1-t)*(thickness_global_data-z_right_inner_global_data) + t * (thickness_global_data-z_right_outer_global_data)))));


    if (std::isnan(z))
    {
        z = 0;
    }


    return D3vector ( x,y,z);

}

D3vector transform_right_lens_inner_quad ( double t1, double t2, double* pointer_global_data) {
    double R_global_data = pointer_global_data[0];
    double r_global_data              = pointer_global_data[1];
    double curvatureLeft_global_data  = pointer_global_data[2];
    double curvatureRight_global_data = pointer_global_data[3];
    double thickness_global_data          = pointer_global_data[4];
    double offsetZ_global_data        = pointer_global_data[5];
    double z_left_inner_global_data   = pointer_global_data[6];
    double z_left_outer_global_data   = pointer_global_data[7];
    double z_right_inner_global_data  = pointer_global_data[8];
    double z_right_outer_global_data  = pointer_global_data[9];
    double temp = t1;
    t1 = t2;
    t2 = temp;
    double t = t2;
    t2 = (1-t2);

    double sign = (curvatureRight_global_data > 0) ? 1 : ((curvatureRight_global_data < 0) ? -1 : 0)  ;


    double actualX = + r_global_data * t1 - r_global_data * t2;
    double actualY = - r_global_data * t1 - r_global_data * t2 + r_global_data;

    actualY = -actualY;
    actualX = -actualX;

    double radiusSquared = (actualY*actualY+actualX*actualX);

    double z = offsetZ_global_data
            +sign*(( sqrt(pow(curvatureRight_global_data,2)- radiusSquared) //pow(r_global_data,2) * (t*t + (1-t)*(1-t)))
                     +sign*(- curvatureRight_global_data + (thickness_global_data-z_right_inner_global_data)) ));
    if (( t1 == 0 && t2 == 0)||( t1 == 1 && t2 == 0)||( t1 == 0 && t2 == 1)||( t1 == 1 && t2 == 1))
    {
        //cout << "x " << actualX << " y " << actualY <<endl;
//        if ( fabs(z) > 1e-8)
//        {
//            cout << " stop  "<<endl;
//        }
    }


    //return D3vector (0 , 0 , fabs((t1-0)*(t2-0)*(t1-1)*(t2-1)) * 10);
        return D3vector (0 , 0 , z);
}

D3vector transform_left_lens_inner_quad ( double t1, double t2, double* pointer_global_data) {
    double R_global_data              = pointer_global_data[0];
    double r_global_data              = pointer_global_data[1];
    double curvatureLeft_global_data  = pointer_global_data[2];
    double curvatureRight_global_data = pointer_global_data[3];
    double thickness_global_data      = pointer_global_data[4];
    double offsetZ_global_data        = pointer_global_data[5];
    double z_left_inner_global_data   = pointer_global_data[6];
    double z_left_outer_global_data   = pointer_global_data[7];
    double z_right_inner_global_data  = pointer_global_data[8];
    double z_right_outer_global_data  = pointer_global_data[9];
    double temp = t1;
    t1 = t2;
    t2 = temp;
    double t = t2;
    t2 = (1-t2);



    double sign = (curvatureLeft_global_data > 0) ? 1 : ((curvatureLeft_global_data < 0) ? -1 : 0)  ;


    double actualX = + r_global_data * t1 - r_global_data * t2;
    double actualY = - r_global_data * t1 - r_global_data * t2 + r_global_data;

    actualY = -actualY;
    actualX = -actualX;

    double actualRadius = (actualY*actualY+actualX*actualX);

    double z = 0;


    z = offsetZ_global_data+sign*(-1*( sqrt(pow(curvatureLeft_global_data,2)-actualRadius) - sign* (curvatureLeft_global_data - z_left_inner_global_data)) );

    if (( t1 == 0 && t2 == 0)||( t1 == 1 && t2 == 0)||( t1 == 0 && t2 == 1)||( t1 == 1 && t2 == 1))
    {
        //cout << "x " << actualX << " y " << actualY <<endl;
//        if ( fabs(z) > 1e-8)
//        {
//            cout << " stop  "<<endl;
//        }
    }
        return D3vector (0 , 0 , z);

}

D3vector transform_left_lens_diag_NW_quad ( double t1, double t2, double* pointer_global_data) {
    double R_global_data = pointer_global_data[0];
    double r_global_data              = pointer_global_data[1];
    double curvatureLeft_global_data  = pointer_global_data[2];
    double curvatureRight_global_data = pointer_global_data[3];
    double thickness_global_data      = pointer_global_data[4];
    double offsetZ_global_data        = pointer_global_data[5];
    double z_left_inner_global_data   = pointer_global_data[6];
    double z_left_outer_global_data   = pointer_global_data[7];
    double z_right_inner_global_data  = pointer_global_data[8];
    double z_right_outer_global_data  = pointer_global_data[9];
    double temp = t1;
    t1 = t2;
    t2 = temp;
    double t = t2;
   // t2 = (1-t2);
    double x = -R_global_data*( sin ( t1*0.5*M_PI )-t1);
    x = x*t2;
    double y = R_global_data*(cos ( t1*0.5*M_PI )- ( 1-t1 )) ;
    y = y*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;

    actualX = -actualX;
    double radiusSquared = ((x+actualX)*(x+actualX)+(y+actualY)*(y+actualY));

    double sign = (curvatureLeft_global_data > 0) ? 1 : ((curvatureLeft_global_data < 0) ? -1 : 0)  ;

    double z =offsetZ_global_data+sign*(-1*( sqrt(pow(curvatureLeft_global_data,2)-radiusSquared) +sign *(-curvatureLeft_global_data +( (1-t)*z_left_inner_global_data + t * z_left_outer_global_data))));
     if (std::isnan(z))
     {z = 0;}

     return D3vector ( x,y,z);
}





D3vector transform_right_lens_diag_SW_quad ( double t1, double t2, double* pointer_global_data) {
    double R_global_data = pointer_global_data[0];
    double r_global_data              = pointer_global_data[1];
    double curvatureLeft_global_data  = pointer_global_data[2];
    double curvatureRight_global_data = pointer_global_data[3];
    double thickness_global_data      = pointer_global_data[4];
    double offsetZ_global_data        = pointer_global_data[5];
    double z_left_inner_global_data   = pointer_global_data[6];
    double z_left_outer_global_data   = pointer_global_data[7];
    double z_right_inner_global_data  = pointer_global_data[8];
    double z_right_outer_global_data  = pointer_global_data[9];
    double temp = t1;
    t1 = t2;
    t2 = temp;
    double t = t2;
    double x = -R_global_data *( sin ( t1*0.5*M_PI )-t1);
    x = x*t2;
    double y = -R_global_data *(cos ( t1*0.5*M_PI )- ( 1-t1 )) ;
    y = y*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;
    actualY = -actualY;
    actualX = -actualX;
    double sign = (curvatureRight_global_data > 0) ? 1 : ((curvatureRight_global_data < 0) ? -1 : 0)  ;

    double z = offsetZ_global_data+
            sign*(  ( sqrt(pow(curvatureRight_global_data,2)-
                           ((x+actualX)*(x+actualX)+(y+actualY)*(y+actualY)))
                      + sign*(-curvatureRight_global_data +( (1-t)*(thickness_global_data-z_right_inner_global_data) + t * (thickness_global_data-z_right_outer_global_data)))));


    if (std::isnan(z))
    {        z = 0;    }
    return D3vector ( x,y,z);

}

D3vector transform_left_lens_diag_SW_quad ( double t1, double t2, double* pointer_global_data) {
    double R_global_data              = pointer_global_data[0];
    double r_global_data              = pointer_global_data[1];
    double curvatureLeft_global_data  = pointer_global_data[2];
    double curvatureRight_global_data = pointer_global_data[3];
    double thickness_global_data      = pointer_global_data[4];
    double offsetZ_global_data        = pointer_global_data[5];
    double z_left_inner_global_data   = pointer_global_data[6];
    double z_left_outer_global_data   = pointer_global_data[7];
    double z_right_inner_global_data  = pointer_global_data[8];
    double z_right_outer_global_data  = pointer_global_data[9];
    double temp = t1;
    t1 = t2;
    t2 = temp;
    double t = t2;
    double x = -R_global_data *( sin ( t1*0.5*M_PI )-t1);
    x = x*t2;
    double y = -R_global_data *(cos ( t1*0.5*M_PI )- ( 1-t1 )) ;
    y = y*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;
    actualY = -actualY;
    actualX = -actualX;

    double sign = (curvatureLeft_global_data > 0) ? 1 : ((curvatureLeft_global_data < 0) ? -1 : 0)  ;

    double radiusSquared = ((x+actualX)*(x+actualX)+(y+actualY)*(y+actualY));

    double z =offsetZ_global_data+sign*(-1*( sqrt(pow(curvatureLeft_global_data,2)-radiusSquared) +sign *(-curvatureLeft_global_data +( (1-t)*z_left_inner_global_data + t * z_left_outer_global_data))));
     if (std::isnan(z))
     {z = 0;}

     return D3vector ( x,y,z);
}


D3vector transform_left_lens_diag_SE_quad ( double t1, double t2, double* pointer_global_data) {
    double R_global_data              = pointer_global_data[0];
    double r_global_data              = pointer_global_data[1];
    double curvatureLeft_global_data  = pointer_global_data[2];
    double curvatureRight_global_data = pointer_global_data[3];
    double thickness_global_data      = pointer_global_data[4];
    double offsetZ_global_data        = pointer_global_data[5];
    double z_left_inner_global_data   = pointer_global_data[6];
    double z_left_outer_global_data   = pointer_global_data[7];
    double z_right_inner_global_data  = pointer_global_data[8];
    double z_right_outer_global_data  = pointer_global_data[9];
    double temp = t1;
    t1 = t2;
    t2 = temp;
    double t = t2;
    double x =R_global_data * ( sin ( t1*0.5*M_PI )-t1);
    x = x*t2;
    double y = -R_global_data *(cos ( t1*0.5*M_PI )- ( 1-t1 )) ;
    y = y*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;
    actualY = -actualY;

    double sign = (curvatureLeft_global_data > 0) ? 1 : ((curvatureLeft_global_data < 0) ? -1 : 0)  ;
    double radiusSquared = ((x+actualX)*(x+actualX)+(y+actualY)*(y+actualY));

    //std::cout << "radius inner " << sqrt(radiusSquared ) << std::endl;
    double z =offsetZ_global_data+sign*(-1*( sqrt(pow(curvatureLeft_global_data,2)-radiusSquared) +sign *(-curvatureLeft_global_data +( (1-t)*z_left_inner_global_data + t * z_left_outer_global_data))));
     if (std::isnan(z))
     {z = 0;}

     return D3vector ( x,y,z);
}

D3vector transform_right_lens_diag_SE_quad ( double t1, double t2, double* pointer_global_data) {
    double R_global_data              = pointer_global_data[0];
    double r_global_data              = pointer_global_data[1];
    double curvatureLeft_global_data  = pointer_global_data[2];
    double curvatureRight_global_data = pointer_global_data[3];
    double thickness_global_data      = pointer_global_data[4];
    double offsetZ_global_data        = pointer_global_data[5];
    double z_left_inner_global_data   = pointer_global_data[6];
    double z_left_outer_global_data   = pointer_global_data[7];
    double z_right_inner_global_data  = pointer_global_data[8];
    double z_right_outer_global_data  = pointer_global_data[9];
    double temp = t1;
    t1 = t2;
    t2 = temp;
    double t = t2;
    double x =R_global_data * ( sin ( t1*0.5*M_PI )-t1);
    x = x*t2;
    double y = -R_global_data *(cos ( t1*0.5*M_PI )- ( 1-t1 )) ;
    y = y*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;
    actualY = -actualY;
    double sign = (curvatureRight_global_data > 0) ? 1 : ((curvatureRight_global_data < 0) ? -1 : 0)  ;

    double z = offsetZ_global_data+
            sign*(  ( sqrt(pow(curvatureRight_global_data,2)-
                           ((x+actualX)*(x+actualX)+(y+actualY)*(y+actualY)))
                      + sign*(-curvatureRight_global_data +( (1-t)*(thickness_global_data-z_right_inner_global_data) + t * (thickness_global_data-z_right_outer_global_data)))));


    if (std::isnan(z))
    {
    z = 0;}
    return D3vector ( x,y,z);

}

D3vector transform_left_lens_diag_NE_quad ( double t1, double t2, double* pointer_global_data) {
    double R_global_data              = pointer_global_data[0];
    double r_global_data              = pointer_global_data[1];
    double curvatureLeft_global_data  = pointer_global_data[2];
    double curvatureRight_global_data = pointer_global_data[3];
    double thickness_global_data      = pointer_global_data[4];
    double offsetZ_global_data        = pointer_global_data[5];
    double z_left_inner_global_data   = pointer_global_data[6];
    double z_left_outer_global_data   = pointer_global_data[7];
    double z_right_inner_global_data  = pointer_global_data[8];
    double z_right_outer_global_data  = pointer_global_data[9];
    double temp = t1;

    t1 = t2;
    t2 = temp;
    double t = t2;

    double x = R_global_data * ( sin ( t1*0.5*M_PI )-t1);

    x = x*t2;
    double y = R_global_data *  (cos ( t1*0.5*M_PI )- ( 1-t1 )) ;
    y = y*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 sign = (curvatureLeft_global_data > 0) ? 1 : ((curvatureLeft_global_data < 0) ? -1 : 0)  ;

    double radiusSquared = ((x+actualX)*(x+actualX)+(y+actualY)*(y+actualY));


    double z  =offsetZ_global_data+sign*(-1*( sqrt(pow(curvatureLeft_global_data,2)-radiusSquared) +sign *(-curvatureLeft_global_data +( (1-t)*z_left_inner_global_data + t * z_left_outer_global_data))));
    if (std::isnan(z))
    {z = 0;}

    return D3vector ( x,y,z);
}

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[1];
    double curvatureLeft_global_data  = pointer_global_data[2];
    double curvatureRight_global_data = pointer_global_data[3];
    double thickness_global_data      = pointer_global_data[4];
    double offsetZ_global_data        = pointer_global_data[5];
    double z_left_inner_global_data   = pointer_global_data[6];
    double z_left_outer_global_data   = pointer_global_data[7];
    double z_right_inner_global_data  = pointer_global_data[8];
    double z_right_outer_global_data  = pointer_global_data[9];
    bool invertT1T2 = true;
    if (invertT1T2)
    {

        double temp = t1;
        t1 = t2;
        t2 = temp;
    }
    double t = t2;



    double sign = (curvatureRight_global_data > 0) ? 1 : ((curvatureRight_global_data < 0) ? -1 : 0)  ;
    double signLeft = (curvatureLeft_global_data > 0) ? 1 : ((curvatureLeft_global_data < 0) ? -1 : 0)  ;

    double radiusSquared = r_global_data + (R_global_data-r_global_data)* (t2);



    radiusSquared = radiusSquared * radiusSquared;



    double zRight = offsetZ_global_data+
            sign*(  ( sqrt(pow(curvatureRight_global_data,2)-
                          ( radiusSquared))
                      + sign*(-curvatureRight_global_data +( (1-t)*(thickness_global_data-z_right_inner_global_data) + t * (thickness_global_data-z_right_outer_global_data)))));
    double zLeft =   offsetZ_global_data
                    +signLeft*(-1*( sqrt(pow(curvatureLeft_global_data,2)-
                            radiusSquared)
                        +signLeft *(-curvatureLeft_global_data +( (1-t)*z_left_inner_global_data + t * z_left_outer_global_data))));



   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 (  fabs(z) > 1e-8)
//        {
//            cout << "stop  "<<endl;
//        }
    }
        return D3vector (   0 ,0 ,z);

}

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[1];
    double curvatureLeft_global_data  = pointer_global_data[2];
    double curvatureRight_global_data = pointer_global_data[3];
    double thickness_global_data      = pointer_global_data[4];
    double offsetZ_global_data        = pointer_global_data[5];
    double z_left_inner_global_data   = pointer_global_data[6];
    double z_left_outer_global_data   = pointer_global_data[7];
    double z_right_inner_global_data  = pointer_global_data[8];
    double z_right_outer_global_data  = pointer_global_data[9];

    double z_right_edge_global_data    = pointer_global_data[10];
    double z_left_edge_global_data     = pointer_global_data[11];
    double R_outer_global_data        = pointer_global_data[12];




    bool invertT1T2 = true;
    if (invertT1T2)
    {

        double temp = t1;
        t1 = t2;
        t2 = temp;
    }
    double t = t2;



    double sign = (curvatureRight_global_data > 0) ? 1 : ((curvatureRight_global_data < 0) ? -1 : 0)  ;
    double signLeft = (curvatureLeft_global_data > 0) ? 1 : ((curvatureLeft_global_data < 0) ? -1 : 0)  ;

    double radiusSquared = R_global_data + (R_outer_global_data-R_global_data)* (t2);
    radiusSquared = radiusSquared * radiusSquared;

    double zRight = offsetZ_global_data+
            sign*(  ( sqrt(pow(curvatureRight_global_data,2)-
                          ( radiusSquared))
                      + sign*(-curvatureRight_global_data +( (1-t)*(thickness_global_data-z_right_outer_global_data) + t * (thickness_global_data-z_right_edge_global_data)))));

    if (z_right_outer_global_data == z_right_edge_global_data)
    {
        zRight = 0;
    }
    double zLeft =   offsetZ_global_data
                    +signLeft*(-1*( sqrt(pow(curvatureLeft_global_data,2)-
                            radiusSquared)
                        +signLeft *(-curvatureLeft_global_data +( (1-t)*z_left_outer_global_data + t * z_left_edge_global_data))));

    if (z_left_outer_global_data == z_left_edge_global_data)
    {
        zLeft = 0;
    }


   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 (  fabs(z) > 1e-8)
//        {
//            cout << "stop  "<<endl;
//        }
    }
        return D3vector (   0 ,0 ,z);

}


D3vector transform_inner_faces_NE_quad( double t1, double t2, double* pointer_global_data) {
    double R_global_data              = pointer_global_data[0];
    double r_global_data              = pointer_global_data[1];
    double curvatureLeft_global_data  = pointer_global_data[2];
    double curvatureRight_global_data = pointer_global_data[3];
    double thickness_global_data      = pointer_global_data[4];
    double offsetZ_global_data        = pointer_global_data[5];
    double z_left_inner_global_data   = pointer_global_data[6];
    double z_left_outer_global_data   = pointer_global_data[7];
    double z_right_inner_global_data  = pointer_global_data[8];
    double z_right_outer_global_data  = pointer_global_data[9];
    bool invertT1T2 = false;
    if (invertT1T2)
    {

        double temp = t1;
        t1 = t2;
        t2 = temp;
    }

    double sign = (curvatureRight_global_data > 0) ? 1 : ((curvatureRight_global_data < 0) ? -1 : 0)  ;
    double signLeft = (curvatureLeft_global_data > 0) ? 1 : ((curvatureLeft_global_data < 0) ? -1 : 0)  ;

    double x = r_global_data* (t1-1);
    double y = r_global_data* (t1);
    double radiusSquared = x*x+y*y;

    double zRight = offsetZ_global_data
                +sign*(( sqrt(pow(curvatureRight_global_data,2)- radiusSquared) //pow(r_global_data,2) * (t*t + (1-t)*(1-t)))
                         +sign*(- curvatureRight_global_data + (thickness_global_data-z_right_inner_global_data)) ));
    double zLeft = offsetZ_global_data
                +signLeft*(-1*( sqrt(pow(curvatureLeft_global_data,2)- radiusSquared)
                            - signLeft* (curvatureLeft_global_data - z_left_inner_global_data)) );

    //works for one side only:
    double z = zRight * (t2) + zLeft * (1-t2);

    if (( t1 == 0 && t2 == 0)||( t1 == 1 && t2 == 0)||( t1 == 0 && t2 == 1)||( t1 == 1 && t2 == 1))
    {
//        if (  fabs(z) > 1e-8)
//        {
//            cout << "stop  "<<endl;
//        }
    }
        return D3vector (   0 ,0 ,z);

}

D3vector transform_right_lens_diag_NE_quad ( double t1, double t2, double* pointer_global_data) {
    double R_global_data              = pointer_global_data[0];
    double r_global_data              = pointer_global_data[1];
    double curvatureLeft_global_data  = pointer_global_data[2];
    double curvatureRight_global_data = pointer_global_data[3];
    double thickness_global_data      = pointer_global_data[4];
    double offsetZ_global_data        = pointer_global_data[5];
    double z_left_inner_global_data   = pointer_global_data[6];
    double z_left_outer_global_data   = pointer_global_data[7];
    double z_right_inner_global_data  = pointer_global_data[8];
    double z_right_outer_global_data  = pointer_global_data[9];
    double temp = t1;
    t1 = t2;
    t2 = temp;
    double t = t2;

    double x = R_global_data * ( sin ( t1*0.5*M_PI )-t1);

    x = x*t2;
    double y = R_global_data *  (cos ( t1*0.5*M_PI )- ( 1-t1 )) ;
    y = y*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 sign = (curvatureRight_global_data > 0) ? 1 : ((curvatureRight_global_data < 0) ? -1 : 0)  ;

    double radiusSquared = ((x+actualX)*(x+actualX)+(y+actualY)*(y+actualY));



    double z = offsetZ_global_data+
            sign*(  ( sqrt(pow(curvatureRight_global_data,2)-
                           radiusSquared)
                      + 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))
    {   z = 0;}
//    std::cout << "transform_right_lens_diag_NE_quad ";
//    D3vector ( x,y,z).Print();
//    std::cout << std::endl;
    return D3vector ( x,y,z);

}


D3vector transform_right_lens_diag_NW_quad_cut ( double t1, double t2, double* pointer_global_data) {

    double R_global_data              = pointer_global_data[0];
    double r_global_data              = pointer_global_data[1];
    double curvatureLeft_global_data  = pointer_global_data[2];
    double curvatureRight_global_data = pointer_global_data[3];
    double thickness_global_data      = pointer_global_data[4];
    double offsetZ_global_data        = pointer_global_data[5];
    double z_left_inner_global_data   = pointer_global_data[6];
    double z_left_outer_global_data   = pointer_global_data[7];
    double z_right_inner_global_data  = pointer_global_data[8];
    double z_right_outer_global_data  = pointer_global_data[9];

    //shirts the transformation to the cut domain
    r_global_data                     = R_global_data;
    R_global_data                     = pointer_global_data[12];

    z_left_inner_global_data          = z_left_outer_global_data;
    z_right_inner_global_data         = z_right_outer_global_data;

    z_left_outer_global_data          = pointer_global_data[11];
    z_right_outer_global_data         = pointer_global_data[10];


    double temp = t1;
    t1 = t2;
    t2 = temp;
    double t = t2;
   // t2 = (1-t2);
    double x = -R_global_data*( sin ( t1*0.5*M_PI )-t1);
    x = x*t2;
    double y = R_global_data*(cos ( t1*0.5*M_PI )- ( 1-t1 )) ;
    y = y*t2;


    //added due to bended inner part :
    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 )) ;
    x += (1-t2) * xAdd;
    y += (1-t2) * yAdd;

    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;

    actualX = -actualX;
    double radiusSquared = ((x+actualX)*(x+actualX)+(y+actualY)*(y+actualY));
    double sign = (curvatureRight_global_data > 0) ? 1 : ((curvatureRight_global_data < 0) ? -1 : 0)  ;

    double z = offsetZ_global_data+
                sign*(  ( sqrt(pow(curvatureRight_global_data,2)-
                               radiusSquared)
                          + sign*(-curvatureRight_global_data +( (1-t)*(thickness_global_data-z_right_inner_global_data) + t * (thickness_global_data-z_right_outer_global_data)))));

    if (z_right_inner_global_data == z_right_outer_global_data)
    {
        z = 0;
    }

    if (std::isnan(z))
    {
        z = 0;
    }


    return D3vector ( x,y,z);

}

D3vector transform_left_lens_diag_NW_quad_cut ( double t1, double t2, double* pointer_global_data) {
    double R_global_data = pointer_global_data[0];
    double r_global_data              = pointer_global_data[1];
    double curvatureLeft_global_data  = pointer_global_data[2];
    double curvatureRight_global_data = pointer_global_data[3];
    double thickness_global_data      = pointer_global_data[4];
    double offsetZ_global_data        = pointer_global_data[5];
    double z_left_inner_global_data   = pointer_global_data[6];
    double z_left_outer_global_data   = pointer_global_data[7];
    double z_right_inner_global_data  = pointer_global_data[8];
    double z_right_outer_global_data  = pointer_global_data[9];


    //shirts the transformation to the cut domain
    r_global_data                     = R_global_data;
    R_global_data                     = pointer_global_data[12];

    z_left_inner_global_data          = z_left_outer_global_data;
    z_right_inner_global_data         = z_right_outer_global_data;

    z_left_outer_global_data          = pointer_global_data[11];
    z_right_outer_global_data         = pointer_global_data[10];

    double temp = t1;
    t1 = t2;
    t2 = temp;
    double t = t2;
   // t2 = (1-t2);
    double x = -R_global_data*( sin ( t1*0.5*M_PI )-t1);
    x = x*t2;
    double y = R_global_data*(cos ( t1*0.5*M_PI )- ( 1-t1 )) ;
    y = y*t2;

    //added due to bended inner part :
    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 )) ;
    x += (1-t2) * xAdd;
    y += (1-t2) * yAdd;

    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;

    actualX = -actualX;
    double radiusSquared = ((x+actualX)*(x+actualX)+(y+actualY)*(y+actualY));

    double sign = (curvatureLeft_global_data > 0) ? 1 : ((curvatureLeft_global_data < 0) ? -1 : 0)  ;

    double z =offsetZ_global_data+sign*(-1*( sqrt(pow(curvatureLeft_global_data,2)-radiusSquared) +sign *(-curvatureLeft_global_data +( (1-t)*z_left_inner_global_data + t * z_left_outer_global_data))));

    if (z_left_inner_global_data == z_left_outer_global_data)
    {
        z = 0;
    }

    if (std::isnan(z))
     {z = 0;}

     return D3vector ( x,y,z);
}

D3vector transform_right_lens_diag_SW_quad_cut ( double t1, double t2, double* pointer_global_data) {
    double R_global_data = pointer_global_data[0];
    double r_global_data              = pointer_global_data[1];
    double curvatureLeft_global_data  = pointer_global_data[2];
    double curvatureRight_global_data = pointer_global_data[3];
    double thickness_global_data      = pointer_global_data[4];
    double offsetZ_global_data        = pointer_global_data[5];
    double z_left_inner_global_data   = pointer_global_data[6];
    double z_left_outer_global_data   = pointer_global_data[7];
    double z_right_inner_global_data  = pointer_global_data[8];
    double z_right_outer_global_data  = pointer_global_data[9];

    //shirts the transformation to the cut domain
    r_global_data                     = R_global_data;
    R_global_data                     = pointer_global_data[12];

    z_left_inner_global_data          = z_left_outer_global_data;
    z_right_inner_global_data         = z_right_outer_global_data;

    z_left_outer_global_data          = pointer_global_data[11];
    z_right_outer_global_data         = pointer_global_data[10];

    double temp = t1;
    t1 = t2;
    t2 = temp;
    double t = t2;
    double x = -R_global_data *( sin ( t1*0.5*M_PI )-t1);
    x = x*t2;
    double y = -R_global_data *(cos ( t1*0.5*M_PI )- ( 1-t1 )) ;
    y = y*t2;

    //added due to bended inner part :
    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 )) ;
    x += (1-t2) * xAdd;
    y += (1-t2) * yAdd;

    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;
    actualY = -actualY;
    actualX = -actualX;
    double sign = (curvatureRight_global_data > 0) ? 1 : ((curvatureRight_global_data < 0) ? -1 : 0)  ;

    double z = offsetZ_global_data+
            sign*(  ( sqrt(pow(curvatureRight_global_data,2)-
                           ((x+actualX)*(x+actualX)+(y+actualY)*(y+actualY)))
                      + sign*(-curvatureRight_global_data +( (1-t)*(thickness_global_data-z_right_inner_global_data) + t * (thickness_global_data-z_right_outer_global_data)))));

    if (z_right_inner_global_data == z_right_outer_global_data)
    {
        z = 0;
    }
    if (std::isnan(z))
    {        z = 0;    }
    return D3vector ( x,y,z);

}

D3vector transform_left_lens_diag_SW_quad_cut ( double t1, double t2, double* pointer_global_data) {
    double R_global_data              = pointer_global_data[0];
    double r_global_data              = pointer_global_data[1];
    double curvatureLeft_global_data  = pointer_global_data[2];
    double curvatureRight_global_data = pointer_global_data[3];
    double thickness_global_data      = pointer_global_data[4];
    double offsetZ_global_data        = pointer_global_data[5];
    double z_left_inner_global_data   = pointer_global_data[6];
    double z_left_outer_global_data   = pointer_global_data[7];
    double z_right_inner_global_data  = pointer_global_data[8];
    double z_right_outer_global_data  = pointer_global_data[9];

    //shirts the transformation to the cut domain
    r_global_data                     = R_global_data;
    R_global_data                     = pointer_global_data[12];

    z_left_inner_global_data          = z_left_outer_global_data;
    z_right_inner_global_data         = z_right_outer_global_data;

    z_left_outer_global_data          = pointer_global_data[11];
    z_right_outer_global_data         = pointer_global_data[10];

    double temp = t1;
    t1 = t2;
    t2 = temp;
    double t = t2;
    double x = -R_global_data *( sin ( t1*0.5*M_PI )-t1);
    x = x*t2;
    double y = -R_global_data *(cos ( t1*0.5*M_PI )- ( 1-t1 )) ;
    y = y*t2;

    //added due to bended inner part :
    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 )) ;
    x += (1-t2) * xAdd;
    y += (1-t2) * yAdd;

    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;
    actualY = -actualY;
    actualX = -actualX;

    double sign = (curvatureLeft_global_data > 0) ? 1 : ((curvatureLeft_global_data < 0) ? -1 : 0)  ;

    double radiusSquared = ((x+actualX)*(x+actualX)+(y+actualY)*(y+actualY));

    double z =offsetZ_global_data+sign*(-1*( sqrt(pow(curvatureLeft_global_data,2)-radiusSquared) +sign *(-curvatureLeft_global_data +( (1-t)*z_left_inner_global_data + t * z_left_outer_global_data))));

    if (z_left_inner_global_data == z_left_outer_global_data)
    {
        z = 0;
    }

    if (std::isnan(z))
     {z = 0;}

     return D3vector ( x,y,z);
}

D3vector transform_left_lens_diag_SE_quad_cut ( double t1, double t2, double* pointer_global_data) {
    double R_global_data              = pointer_global_data[0];
    double r_global_data              = pointer_global_data[1];
    double curvatureLeft_global_data  = pointer_global_data[2];
    double curvatureRight_global_data = pointer_global_data[3];
    double thickness_global_data      = pointer_global_data[4];
    double offsetZ_global_data        = pointer_global_data[5];
    double z_left_inner_global_data   = pointer_global_data[6];
    double z_left_outer_global_data   = pointer_global_data[7];
    double z_right_inner_global_data  = pointer_global_data[8];
    double z_right_outer_global_data  = pointer_global_data[9];

    //shirts the transformation to the cut domain
    r_global_data                     = R_global_data;
    R_global_data                     = pointer_global_data[12];

    z_left_inner_global_data          = z_left_outer_global_data;
    z_right_inner_global_data         = z_right_outer_global_data;

    z_left_outer_global_data          = pointer_global_data[11];
    z_right_outer_global_data         = pointer_global_data[10];

    double temp = t1;
    t1 = t2;
    t2 = temp;
    double t = t2;
    double x =R_global_data * ( sin ( t1*0.5*M_PI )-t1);
    x = x*t2;
    double y = -R_global_data *(cos ( t1*0.5*M_PI )- ( 1-t1 )) ;
    y = y*t2;
    //added due to bended inner part :
    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 )) ;
    x += (1-t2) * xAdd;
    y += (1-t2) * yAdd;

    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;
    actualY = -actualY;

    double sign = (curvatureLeft_global_data > 0) ? 1 : ((curvatureLeft_global_data < 0) ? -1 : 0)  ;
    double radiusSquared = ((x+actualX)*(x+actualX)+(y+actualY)*(y+actualY));
    //std::cout << "radius outer " << sqrt(radiusSquared ) << std::endl;

    double z =offsetZ_global_data+sign*(-1*( sqrt(pow(curvatureLeft_global_data,2)-radiusSquared) +sign *(-curvatureLeft_global_data +( (1-t)*z_left_inner_global_data + t * z_left_outer_global_data))));

    if (z_left_inner_global_data == z_left_outer_global_data)
    {
        z = 0;
    }

    if (std::isnan(z))
     {z = 0;}


     return D3vector ( x,y,z);
}

D3vector transform_right_lens_diag_SE_quad_cut ( double t1, double t2, double* pointer_global_data) {
    double R_global_data              = pointer_global_data[0];
    double r_global_data              = pointer_global_data[1];
    double curvatureLeft_global_data  = pointer_global_data[2];
    double curvatureRight_global_data = pointer_global_data[3];
    double thickness_global_data      = pointer_global_data[4];
    double offsetZ_global_data        = pointer_global_data[5];
    double z_left_inner_global_data   = pointer_global_data[6];
    double z_left_outer_global_data   = pointer_global_data[7];
    double z_right_inner_global_data  = pointer_global_data[8];
    double z_right_outer_global_data  = pointer_global_data[9];

    //shirts the transformation to the cut domain
    r_global_data                     = R_global_data;
    R_global_data                     = pointer_global_data[12];

    z_left_inner_global_data          = z_left_outer_global_data;
    z_right_inner_global_data         = z_right_outer_global_data;

    z_left_outer_global_data          = pointer_global_data[11];
    z_right_outer_global_data         = pointer_global_data[10];

    double temp = t1;
    t1 = t2;
    t2 = temp;
    double t = t2;
    double x =R_global_data * ( sin ( t1*0.5*M_PI )-t1);
    x = x*t2;
    double y = -R_global_data *(cos ( t1*0.5*M_PI )- ( 1-t1 )) ;
    y = y*t2;

    //added due to bended inner part :
    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 )) ;
    x += (1-t2) * xAdd;
    y += (1-t2) * yAdd;

    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;
    actualY = -actualY;
    double sign = (curvatureRight_global_data > 0) ? 1 : ((curvatureRight_global_data < 0) ? -1 : 0)  ;

    double z = offsetZ_global_data+
            sign*(  ( sqrt(pow(curvatureRight_global_data,2)-
                           ((x+actualX)*(x+actualX)+(y+actualY)*(y+actualY)))
                      + sign*(-curvatureRight_global_data +( (1-t)*(thickness_global_data-z_right_inner_global_data) + t * (thickness_global_data-z_right_outer_global_data)))));

    if (z_right_inner_global_data == z_right_outer_global_data)
    {
        z = 0;
    }

    if (std::isnan(z))
    {
    z = 0;}
    return D3vector ( x,y,z);

}

D3vector transform_left_lens_diag_NE_quad_cut ( double t1, double t2, double* pointer_global_data) {
    double R_global_data              = pointer_global_data[0];
    double r_global_data              = pointer_global_data[1];
    double curvatureLeft_global_data  = pointer_global_data[2];
    double curvatureRight_global_data = pointer_global_data[3];
    double thickness_global_data      = pointer_global_data[4];
    double offsetZ_global_data        = pointer_global_data[5];
    double z_left_inner_global_data   = pointer_global_data[6];
    double z_left_outer_global_data   = pointer_global_data[7];
    double z_right_inner_global_data  = pointer_global_data[8];
    double z_right_outer_global_data  = pointer_global_data[9];

    //shirts the transformation to the cut domain
    r_global_data                     = R_global_data;
    R_global_data                     = pointer_global_data[12];

    z_left_inner_global_data          = z_left_outer_global_data;
    z_right_inner_global_data         = z_right_outer_global_data;

    z_left_outer_global_data          = pointer_global_data[11];
    z_right_outer_global_data         = pointer_global_data[10];

    double temp = t1;

    t1 = t2;
    t2 = temp;
    double t = t2;

    double x = R_global_data * ( sin ( t1*0.5*M_PI )-t1);

    x = x*t2;
    double y = R_global_data *  (cos ( t1*0.5*M_PI )- ( 1-t1 )) ;
    y = y*t2;
    double xT = x;
    double yT = y;
    //added due to bended inner part :
    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 )) ;
    x += (1-t2) * xAdd;
    y += (1-t2) * yAdd;


    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 sign = (curvatureLeft_global_data > 0) ? 1 : ((curvatureLeft_global_data < 0) ? -1 : 0)  ;

    double radiusSquared = ((x+actualX)*(x+actualX)+(y+actualY)*(y+actualY));


    double z  =offsetZ_global_data+sign*(-1*( sqrt(pow(curvatureLeft_global_data,2)-radiusSquared) +sign *(-curvatureLeft_global_data +( (1-t)*z_left_inner_global_data + t * 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;
    }

    if (std::isnan(z))
    {z = 0;}

    //return D3vector ( xT,yT,z);
    return D3vector ( x,y,z);
}

D3vector transform_right_lens_diag_NE_quad_cut ( double t1, double t2, double* pointer_global_data) {
    double R_global_data              = pointer_global_data[0];
    double r_global_data              = pointer_global_data[1];
    double curvatureLeft_global_data  = pointer_global_data[2];
    double curvatureRight_global_data = pointer_global_data[3];
    double thickness_global_data      = pointer_global_data[4];
    double offsetZ_global_data        = pointer_global_data[5];
    double z_left_inner_global_data   = pointer_global_data[6];
    double z_left_outer_global_data   = pointer_global_data[7];
    double z_right_inner_global_data  = pointer_global_data[8];
    double z_right_outer_global_data  = pointer_global_data[9];

    //shirts the transformation to the cut domain
    r_global_data                     = R_global_data;
    R_global_data                     = pointer_global_data[12];

    z_left_inner_global_data          = z_left_outer_global_data;
    z_right_inner_global_data         = z_right_outer_global_data;

    z_left_outer_global_data          = pointer_global_data[11];
    z_right_outer_global_data         = pointer_global_data[10];

    double temp = t1;
    t1 = t2;
    t2 = temp;
    double t = t2;

    double x = R_global_data * ( sin ( t1*0.5*M_PI )-t1);

    x = x*t2;
    double y = R_global_data *  (cos ( t1*0.5*M_PI )- ( 1-t1 )) ;
    y = y*t2;

    double xT = x;
    double yT = y;
    //added due to bended inner part :
    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 )) ;
    x += (1-t2) * xAdd;
    y += (1-t2) * yAdd;

    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 sign = (curvatureRight_global_data > 0) ? 1 : ((curvatureRight_global_data < 0) ? -1 : 0)  ;

    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+
            sign*(  ( sqrt(pow(curvatureRight_global_data,2)-
                           radiusSquared)
                      + 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)
    {
        z = 0;
    }

    if (std::isnan(x) || std::isnan(y) || std::isnan(z))
    {   z = 0;}

//    std::cout << "transform_right_lens_diag_NE_quad_cut ";
//    D3vector ( x,y,z).Print();
//    std::cout << std::endl;
//    double actualZ =  offsetZ_global_data+
//            sign*(  ( sqrt(pow(curvatureRight_global_data,2)-
//                           radiusSquared)
//                      + sign*(-curvatureRight_global_data))) ;
//    return D3vector (actualX, actualY,actualZ );
    //return D3vector ( xT,yT,z);
    return D3vector ( x,y,z);
    //return D3vector ( (1-t2) * xAdd,(1-t2) * yAdd,z);
    return D3vector ( 0,0,z);

}





D3vector transform_left_lens_diag ( double t, double* pointer_global_data) {
    double R_global_data              = pointer_global_data[0];
    double r_global_data              = pointer_global_data[1];
    double curvatureLeft_global_data  = pointer_global_data[2];
    double offsetZ_global_data        = pointer_global_data[5];
    double z_left_inner_global_data   = pointer_global_data[6];
    double z_left_outer_global_data   = pointer_global_data[7];

    double sign = (curvatureLeft_global_data > 0) ? 1 : ((curvatureLeft_global_data < 0) ? -1 : 0)  ;

    double z = 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))));
    return D3vector ( 0.0 ,0.0 , z );
}

D3vector transform_left_lens_rectangle ( double t, double* pointer_global_data) {
  double r_global_data              = pointer_global_data[1];
  double curvatureLeft_global_data  = pointer_global_data[2];
  double offsetZ_global_data        =  pointer_global_data[5];
  double z_left_inner_global_data   = pointer_global_data[6];

  //std::cout << "t : " << t << " : " << offsetZ_global_data-1*( sqrt(pow(curvatureLeft_global_data,2)-pow(r_global_data,2) * (t*t + (1-t)*(1-t))) - curvatureLeft_global_data + z_left_inner_global_data)  <<std::endl;
  double sign = (curvatureLeft_global_data > 0) ? 1 : ((curvatureLeft_global_data < 0) ? -1 : 0)  ;
    //std::cout << "transform_left_lens_rectangle : t = " << t << " : " <<offsetZ_global_data+sign*(-1*( sqrt(pow(curvatureLeft_global_data,2)-pow(r_global_data,2) * (t*t + (1-t)*(1-t))) - sign* (curvatureLeft_global_data - z_left_inner_global_data)) ) <<std::endl;
  double z = offsetZ_global_data+sign*(-1*( sqrt(pow(curvatureLeft_global_data,2)-pow(r_global_data,2) * (t*t + (1-t)*(1-t))) - sign* (curvatureLeft_global_data - z_left_inner_global_data)) );
  return D3vector ( 0.0 ,0.0 ,z);
                              //offsetZ_global_data-1*( sqrt(pow(curvatureLeft_global_data,2)-pow(r_global_data,2) * (t*t + (1-t)*(1-t))) - curvatureLeft_global_data + z_left_inner_global_data) )
}

D3vector transform_right_lens_diag ( double t, double* pointer_global_data) {
    double R_global_data              = pointer_global_data[0];
    double r_global_data              = pointer_global_data[1];
    double curvatureRight_global_data = pointer_global_data[3];
    double thickness_global_data          = pointer_global_data[4];
    double offsetZ_global_data        =  pointer_global_data[5];
    double z_right_inner_global_data  = pointer_global_data[8];
    double z_right_outer_global_data  = pointer_global_data[9];

    double sign = (curvatureRight_global_data > 0) ? 1 : ((curvatureRight_global_data < 0) ? -1 : 0)  ;
    //std::cout << "transform_right_lens_diag : t = " << t << " : " <<offsetZ_global_data+ sign*( ( sqrt(pow(curvatureRight_global_data,2)-pow(r_global_data * (1-t) + R_global_data * t,2)) + sign*(-curvatureRight_global_data +( (1-t)*(thickness_global_data-z_right_inner_global_data) + t * (thickness_global_data-z_right_outer_global_data))))) <<std::endl;

    double z = offsetZ_global_data+ sign*(  ( sqrt(pow(curvatureRight_global_data,2)-pow(r_global_data * (1-t) + R_global_data * t,2)) + sign*(-curvatureRight_global_data +( (1-t)*(thickness_global_data-z_right_inner_global_data) + t * (thickness_global_data-z_right_outer_global_data)))));

    return D3vector ( 0.0 ,0.0 , z);
}

D3vector transform_right_lens_rectangle ( double t, double* pointer_global_data) {
  double r_global_data              = pointer_global_data[1];
  double curvatureRight_global_data = pointer_global_data[3];
  double thickness_global_data      = pointer_global_data[4];
  double offsetZ_global_data        =  pointer_global_data[5];
  double z_right_inner_global_data  = pointer_global_data[8];

  double sign = (curvatureRight_global_data > 0) ? 1 : ((curvatureRight_global_data < 0) ? -1 : 0)  ;
  //std::cout << "transform_right_lens_rectangle : t = " << t << " : " <<offsetZ_global_data+sign*(( sqrt(pow(curvatureRight_global_data,2)-pow(r_global_data,2) * (t*t + (1-t)*(1-t))) +sign*(- curvatureRight_global_data + (thickness_global_data-z_right_inner_global_data)) )) <<std::endl;

  return D3vector ( 0.0 ,0.0 ,offsetZ_global_data+sign*(( sqrt(pow(curvatureRight_global_data,2)-pow(r_global_data,2) * (t*t + (1-t)*(1-t))) +sign*(- curvatureRight_global_data + (thickness_global_data-z_right_inner_global_data)) )));
}


D3vector transform_lens_NE ( double t, double* pointer_global_data) {
  double R_global_data = pointer_global_data[0];
  return R_global_data*D3vector ( (cos ( t*0.5*M_PI )- ( 1-t )) ,( sin ( t*0.5*M_PI )-t) , 0.0 );
}

D3vector transform_lens_NW ( double t, double* pointer_global_data ) {
  double R_global_data = pointer_global_data[0];
  return R_global_data*D3vector (( -sin ( t*0.5*M_PI ) +t),( cos ( t*0.5*M_PI )- ( 1-t )), 0.0 );
}

D3vector transform_lens_SW ( double t, double* pointer_global_data ) {
  double R_global_data = pointer_global_data[0];
  return R_global_data*D3vector ( -cos ( t*0.5*M_PI ) + ( 1-t ), -sin ( t*0.5*M_PI ) +t, 0.0 );
}

D3vector transform_lens_SE ( double t, double* pointer_global_data ) {
  double R_global_data = pointer_global_data[0];
  return R_global_data*D3vector (( cos ( t*0.5*M_PI )- ( 1-t )),( -sin ( t*0.5*M_PI ) +t), 0.0 );
}

D3vector transform_lens_NE_inner ( double t, double* pointer_global_data) {
  double R_global_data = pointer_global_data[1];
  return R_global_data*D3vector ( cos ( t*0.5*M_PI )- ( 1-t ), sin ( t*0.5*M_PI )-t, 0.0 );
}

D3vector transform_lens_NW_inner ( double t, double* pointer_global_data ) {
  double R_global_data = pointer_global_data[1];
  return R_global_data*D3vector ( -sin ( t*0.5*M_PI ) +t, cos ( t*0.5*M_PI )- ( 1-t ), 0.0 );
}

D3vector transform_lens_SW_inner ( double t, double* pointer_global_data ) {
  double R_global_data = pointer_global_data[1];
  return R_global_data*D3vector ( -cos ( t*0.5*M_PI ) + ( 1-t ), -sin ( t*0.5*M_PI ) +t, 0.0 );
}

D3vector transform_lens_SE_inner ( double t, double* pointer_global_data ) {
  double R_global_data = pointer_global_data[1];
  return R_global_data*D3vector ( cos ( t*0.5*M_PI )- ( 1-t ), -sin ( t*0.5*M_PI ) +t, 0.0 );
}

// Fuer DEformation der 4 inneren Kanten t=0,...,1
D3vector transform_lens_NE_square ( double t, double* pointer_global_data) {
  double r = pointer_global_data[1];
  double rho = pointer_global_data[2];
  //return D3vector(0.0,0.0,0.0); //Keine Aenderung
  return D3vector (rho*(cos(M_PI/6.0*(1.0+t))-0.5)-r*(1.0-t), rho*(sin(M_PI/6.0*(1.0+t))-0.5)-r*t,0.0);
}

D3vector transform_lens_NW_square ( double t, double* pointer_global_data ) {
  double r = pointer_global_data[1];
  double rho = pointer_global_data[2];
  //return D3vector(0.0,0.0,0.0); //Keine Aenderung
  return D3vector (rho*(-sin(M_PI/6.0*(1.0+t))+0.5)+r*t, rho*(cos(M_PI/6.0*(1.0+t))-0.5)-r*(1.0-t),0.0);
}

D3vector transform_lens_SW_square ( double t, double* pointer_global_data ) {
  double r = pointer_global_data[1];
  double rho = pointer_global_data[2];
  //return D3vector(0.0,0.0,0.0); //Keine Aenderung
  return D3vector (-rho*(cos(M_PI/6.0*(1.0+t))-0.5)+r*(1.0-t), rho*(-sin(M_PI/6.0*(1.0+t))+0.5)+r*t,0.0);
}

D3vector transform_lens_SE_square ( double t, double* pointer_global_data ) {
  double r = pointer_global_data[1];
  double rho = pointer_global_data[2];
  //return D3vector(0.0,0.0,0.0); //Keine Aenderung
  return D3vector (rho*(cos(M_PI/6.0*(1.0+t))-0.5)-r*(1.0-t), rho*(-sin(M_PI/6.0*(1.0+t))+0.5)+r*t,0.0);
}


Lens_Geometry::Lens_Geometry (double Radius, double thickness, double curvatureLeft, double curvatureRight, double offsetX, double offsetY, double offsetZ, bool inner_grid_arched, double radius)
{
    if (radius == 0)
    {
        radius = Radius / 2.0;
    }
    assert(Radius > radius);
    size_pointer_global_data = 10;
    if (fabs(curvatureLeft) >1e10)
    {
        curvatureLeft = 1e10;
    }
    if (fabs(curvatureRight) >1e10)
    {
        curvatureRight = 1e10;
    }

    double inversFocal = (1.5 - 1.0) * ( (1.0 / curvatureLeft) - ( 1.0 / -curvatureRight ) + ( 1.5 - 1.0) * thickness / ( 1.5 * -curvatureRight * curvatureLeft ));

   // cout << "Focal length of Lens is " << 1.0 / inversFocal << " (assume n = 1.5 ) ." << endl;
    pointer_global_data = new double[size_pointer_global_data];
    pointer_global_data[0] = Radius;
    pointer_global_data[1] = radius;
    pointer_global_data[2] = curvatureLeft;
    pointer_global_data[3] = curvatureRight;
    pointer_global_data[4] = thickness;
    pointer_global_data[5] = offsetZ;

    assert ((thickness > 0.0) &&  "Thickness must be positive.");
    assert ((Radius > 0.0) &&  "Radius must be positive.");
    assert ((radius > 0.0) &&  "radius must be positive.");



    constructionParameters.resize(8);
    constructionParameters[0] = Radius;
    constructionParameters[1] = radius;
    constructionParameters[2] = thickness;
    constructionParameters[3] = curvatureLeft;
    constructionParameters[4] = curvatureRight;
    constructionParameters[5] = offsetX;
    constructionParameters[6] = offsetY;
    constructionParameters[7] = offsetZ;
    constructionBoolArched = inner_grid_arched;
    // maybe adjust, if lens is negative (most thick part not in the centre but at the edges!)
    relativeCoordVector.push_back(0.0);
    relativeCoordVector.push_back(thickness);
    typFuerSlice = lensUgTyp;


    int num_blocks_z_dir = 1;
    assert(num_blocks_z_dir>0);
    int num_blocks_total = 5 * num_blocks_z_dir;
    int num_points_per_face = 8;
    int numPoints = num_points_per_face*(num_blocks_z_dir+1) ;

    double z_left_inner  = -sqrt( pow(curvatureLeft ,2) - radius * radius ) + curvatureLeft  + offsetZ;
    double z_left_outer  = -sqrt( pow(curvatureLeft ,2) - Radius * Radius ) + curvatureLeft  + offsetZ;
    double z_right_inner = +sqrt( pow(curvatureRight,2) - radius * radius ) - curvatureRight + offsetZ + thickness;
    double z_right_outer = +sqrt( pow(curvatureRight,2) - Radius * Radius ) - curvatureRight + offsetZ + thickness;

    if (curvatureLeft< 0.0)
    {
        z_left_inner  = +sqrt( pow(curvatureLeft ,2) - radius * radius ) + curvatureLeft  + offsetZ;
        z_left_outer  = +sqrt( pow(curvatureLeft ,2) - Radius * Radius ) + curvatureLeft  + offsetZ;
    }
    if (curvatureRight< 0.0)
    {
        z_right_inner = -sqrt( pow(curvatureRight,2) - radius * radius ) - curvatureRight + offsetZ + thickness;
        z_right_outer = -sqrt( pow(curvatureRight,2) - Radius * Radius ) - curvatureRight + offsetZ + thickness;
    }

    if (curvatureLeft == 0.0)
    {
        z_left_inner  = 0.0  + offsetZ;
        z_left_outer  = 0.0  + offsetZ;
    }
    if (curvatureRight == 0.0)
    {
        z_right_inner = 0.0 + offsetZ + thickness;
        z_right_outer = 0.0 + offsetZ + thickness;
    }




//    cout << z_right_outer - z_left_outer << endl;
//    cout << thickness << endl;
//    double diff = z_right_outer - z_left_outer - thickness ;
//    cout <<fabs(z_right_outer - z_left_outer) - thickness << endl;
//    cout << ((fabs(z_right_outer - z_left_outer) - thickness) < 0.0) << endl;
//    double epsilon = 1e-5;
//    double diff2 = fabs(z_right_outer - z_left_outer) - thickness-epsilon;
    if (curvatureRight < 0 && curvatureLeft < 0)
    {
        assert (((fabs(z_right_outer - z_left_outer) - thickness) < 0.0) &&  "Lens has negative thickness at edge! Make lens more thick or increase radius of curvature.");

    }
    if (curvatureRight > 0 && curvatureLeft > 0)
    {
        assert ((((z_right_outer - z_left_outer) ) > 0.0) &&  "Lens has negative thickness at edge! Make lens more thick or increase radius of curvature.");

    }


    pointer_global_data[6] = z_left_inner;
    pointer_global_data[7] = z_left_outer;
    pointer_global_data[8] = z_right_inner;
    pointer_global_data[9] = z_right_outer;


    Set_number_points ( numPoints );
    Set_number_hexahedra ( num_blocks_total );


    Set_coordinate_point( 0, D3vector ( offsetX+radius,offsetY,z_left_inner ));
    Set_coordinate_point( 1, D3vector ( offsetX+Radius,offsetY,z_left_outer ));

    Set_coordinate_point( 2, D3vector ( offsetX,offsetY+radius,z_left_inner ));
    Set_coordinate_point( 3, D3vector ( offsetX,offsetY+Radius,z_left_outer ));

    Set_coordinate_point( 4, D3vector ( offsetX-radius,offsetY,z_left_inner ));
    Set_coordinate_point( 5, D3vector ( offsetX-Radius,offsetY,z_left_outer ));

    Set_coordinate_point( 6, D3vector ( offsetX,offsetY-radius,z_left_inner ));
    Set_coordinate_point( 7, D3vector ( offsetX,offsetY-Radius,z_left_outer ));

    Set_coordinate_point( 8, D3vector ( offsetX+radius,offsetY,z_right_inner ));
    Set_coordinate_point( 9, D3vector ( offsetX+Radius,offsetY,z_right_outer ));

    Set_coordinate_point( 10, D3vector ( offsetX,offsetY+radius,z_right_inner ));
    Set_coordinate_point( 11, D3vector ( offsetX,offsetY+Radius,z_right_outer ));

    Set_coordinate_point( 12, D3vector ( offsetX-radius,offsetY,z_right_inner ));
    Set_coordinate_point( 13, D3vector ( offsetX-Radius,offsetY,z_right_outer ));

    Set_coordinate_point( 14, D3vector ( offsetX,offsetY-radius,z_right_inner ));
    Set_coordinate_point( 15, D3vector ( offsetX,offsetY-Radius,z_right_outer ));

    int num_blocks = 1;
    for ( int i = 0; i < num_blocks; ++i ) {
      Set_hexahedron ( i, i*num_points_per_face + 0, i*num_points_per_face + 2,
                          i*num_points_per_face + 6, i*num_points_per_face + 4,
                          ( i+1 ) * num_points_per_face + 0, ( i+1 ) * num_points_per_face +2,
                          ( i+1 ) * num_points_per_face + 6, ( i+1 ) * num_points_per_face +4 );
    }
    for ( int j = 0; j < num_points_per_face; j=j+2) {
      for ( int i = 0; i < num_blocks; ++i) {
        int myIndex = (j/2+1)*num_blocks +i;
        Set_hexahedron ( myIndex,
                         i*num_points_per_face     +  j+0,
                         i*num_points_per_face     +  j+1,
                         i*num_points_per_face     +((j+2) % num_points_per_face),
                         i*num_points_per_face     +((j+3) % num_points_per_face),
                         (i+1)*num_points_per_face +  j+0,
                         (i+1)*num_points_per_face +  j+1,
                         (i+1)*num_points_per_face +((j+2) % num_points_per_face),
                         (i+1)*num_points_per_face +((j+3) % num_points_per_face)
                       );
     }
    }
  construction_hexahedron_points_done();
  int shift;
  for ( int i = 0; i <= num_blocks_z_dir; ++i ) {
    shift = i * num_points_per_face;
    Set_transformation_edge ( 1+shift, 3+shift, transform_lens_NE );
    Set_transformation_edge ( 3+shift, 5+shift, transform_lens_NW );
    Set_transformation_edge ( 5+shift, 7+shift, transform_lens_SW );
    Set_transformation_edge ( 1+shift, 7+shift, transform_lens_SE );

    if (inner_grid_arched == true) {
      Set_transformation_edge ( 0+shift, 2+shift, transform_lens_NE_square );
      Set_transformation_edge ( 2+shift, 4+shift, transform_lens_NW_square );
      Set_transformation_edge ( 4+shift, 6+shift, transform_lens_SW_square );
      Set_transformation_edge ( 0+shift, 6+shift, transform_lens_SE_square );
    }
  }


  Set_transformation_edge(0,2,transform_left_lens_rectangle);
  Set_transformation_edge(2,4,transform_left_lens_rectangle);
  Set_transformation_edge(4,6,transform_left_lens_rectangle);
  Set_transformation_edge(0,6,transform_left_lens_rectangle);

  Set_transformation_edge(0,1,transform_left_lens_diag);
  Set_transformation_edge(2,3,transform_left_lens_diag);
  Set_transformation_edge(4,5,transform_left_lens_diag);
  Set_transformation_edge(6,7,transform_left_lens_diag);
  Set_transformation_edge(0+num_points_per_face,2+num_points_per_face,transform_right_lens_rectangle);
  Set_transformation_edge(2+num_points_per_face,4+num_points_per_face,transform_right_lens_rectangle);
  Set_transformation_edge(4+num_points_per_face,6+num_points_per_face,transform_right_lens_rectangle);
  Set_transformation_edge(0+num_points_per_face,6+num_points_per_face,transform_right_lens_rectangle);
  Set_transformation_edge(0+num_points_per_face,1+num_points_per_face,transform_right_lens_diag);
  Set_transformation_edge(2+num_points_per_face,3+num_points_per_face,transform_right_lens_diag);
  Set_transformation_edge(4+num_points_per_face,5+num_points_per_face,transform_right_lens_diag);
  Set_transformation_edge(6+num_points_per_face,7+num_points_per_face,transform_right_lens_diag);

  construction_done();
}

Lens_Geometry_Quad::Lens_Geometry_Quad(double Radius, double thickness, double curvatureLeft, double curvatureRight, double offsetX, double offsetY, double offsetZ, bool inner_grid_arched, double radius, double cut_edge_from_left, double cut_edge_from_right)
{



    if (radius == 0)
    {
        radius = Radius / 2.0;
    }
    assert(Radius > radius);
    size_pointer_global_data = 12;
    if (fabs(curvatureLeft) >1e10)
    {
        curvatureLeft = 1e10;
    }
    if (fabs(curvatureRight) >1e10)
    {
        curvatureRight = 1e10;
    }

    double inversFocal = (1.5 - 1.0) * ( (1.0 / curvatureLeft) - ( 1.0 / -curvatureRight ) + ( 1.5 - 1.0) * thickness / ( 1.5 * -curvatureRight * curvatureLeft ));

    focalLength_ = 1.0 / inversFocal;
  //  cout << "Focal length of Lens is " << focalLength_ << " (assume n = 1.5 ) ." << endl;


    thickness_ = thickness;
    offsetX_ = offsetX;
    offsetY_ = offsetY;
    offsetZ_ = offsetZ;

    pointer_global_data = new double[size_pointer_global_data];
    pointer_global_data[0] = Radius;
    pointer_global_data[1] = radius;
    pointer_global_data[2] = curvatureLeft;
    pointer_global_data[3] = curvatureRight;
    pointer_global_data[4] = thickness;
    pointer_global_data[5] = offsetZ;

    assert ((thickness > 0.0) &&  "Thickness must be positive.");
    assert ((Radius > 0.0) &&  "Radius must be positive.");
    assert ((radius > 0.0) &&  "radius must be positive.");



    constructionParameters.resize(12);
    constructionParameters[0] = Radius;
    constructionParameters[1] = radius;
    constructionParameters[2] = thickness;
    constructionParameters[3] = curvatureLeft;
    constructionParameters[4] = curvatureRight;
    constructionParameters[5] = offsetX;
    constructionParameters[6] = offsetY;
    constructionParameters[7] = offsetZ;
    constructionBoolArched = inner_grid_arched;
    // maybe adjust, if lens is negative (most thick part not in the centre but at the edges!)
    relativeCoordVector.push_back(0.0);
    relativeCoordVector.push_back(thickness);
    typFuerSlice = lensUgTyp;

    int num_blocks_z_dir = 1;
    assert(num_blocks_z_dir>0);
    int num_blocks_total = 5 * num_blocks_z_dir;
    int num_points_per_face = 8;
    int numPoints = num_points_per_face*(num_blocks_z_dir+1) ;

    double z_left_inner  = -sqrt( pow(curvatureLeft ,2) - radius * radius ) + curvatureLeft  + offsetZ;
    double z_left_outer  = -sqrt( pow(curvatureLeft ,2) - Radius * Radius ) + curvatureLeft  + offsetZ;
    double z_right_inner = +sqrt( pow(curvatureRight,2) - radius * radius ) - curvatureRight + offsetZ + thickness;
    double z_right_outer = +sqrt( pow(curvatureRight,2) - Radius * Radius ) - curvatureRight + offsetZ + thickness;

    if (curvatureLeft< 0.0)
    {
        z_left_inner  = +sqrt( pow(curvatureLeft ,2) - radius * radius ) + curvatureLeft  + offsetZ;
        z_left_outer  = +sqrt( pow(curvatureLeft ,2) - Radius * Radius ) + curvatureLeft  + offsetZ;
    }
    if (curvatureRight< 0.0)
    {
        z_right_inner = -sqrt( pow(curvatureRight,2) - radius * radius ) - curvatureRight + offsetZ + thickness;
        z_right_outer = -sqrt( pow(curvatureRight,2) - Radius * Radius ) - curvatureRight + offsetZ + thickness;
    }

    if (curvatureLeft == 0.0)
    {
        z_left_inner  = 0.0  + offsetZ;
        z_left_outer  = 0.0  + offsetZ;
    }
    if (curvatureRight == 0.0)
    {
        z_right_inner = 0.0 + offsetZ + thickness;
        z_right_outer = 0.0 + offsetZ + thickness;
    }

    constructionParameters[8] = std::min(offsetZ,z_left_outer);
    constructionParameters[9] = std::max(offsetZ + thickness,z_right_outer);
    constructionParameters[10] = 0;
    constructionParameters[11] = 0;

    boundaryPoints.push_back(D3vector{-Radius-1,-Radius-1,constructionParameters[8]-1});
    boundaryPoints.push_back(D3vector{+Radius+1,+Radius+1,constructionParameters[9]+1});




//    if (curvatureRight < 0 && curvatureLeft < 0)
//    {
//        assert (((fabs(z_right_outer - z_left_outer) - thickness) < 0.0) &&  "Lens has negative thickness at edge! Make lens more thick or increase radius of curvature.");

//    }
    if (curvatureRight > 0 && curvatureLeft > 0)
    {
        assert ((((z_right_outer - z_left_outer) ) > 0.0) &&  "Lens has negative thickness at edge! Make lens more thick or increase radius of curvature.");

    }


    pointer_global_data[6] = z_left_inner;
    pointer_global_data[7] = z_left_outer;
    pointer_global_data[8] = z_right_inner;
    pointer_global_data[9] = z_right_outer;

    std::cout << "adding cutted edges : \n";
    pointer_global_data[10] = z_left_outer + cut_edge_from_left;
    pointer_global_data[11] = z_right_outer - cut_edge_from_right;


//    cout << "fixed lens param., channge again " << endl;
//    radius = 1;
//    Radius = 2;
//    z_left_inner = z_left_outer = 0;
//    z_right_inner = z_right_outer = 1;



    Set_number_points ( numPoints );
    Set_number_hexahedra ( num_blocks_total );


    Set_coordinate_point( 0, D3vector ( offsetX+radius,offsetY,z_left_inner ));
    Set_coordinate_point( 1, D3vector ( offsetX+Radius,offsetY,z_left_outer ));

    Set_coordinate_point( 2, D3vector ( offsetX       ,offsetY+radius,z_left_inner ));
    Set_coordinate_point( 3, D3vector ( offsetX       ,offsetY+Radius,z_left_outer ));

    Set_coordinate_point( 4, D3vector ( offsetX-radius,offsetY,z_left_inner ));
    Set_coordinate_point( 5, D3vector ( offsetX-Radius,offsetY,z_left_outer ));

    Set_coordinate_point( 6, D3vector ( offsetX,offsetY-radius,z_left_inner ));
    Set_coordinate_point( 7, D3vector ( offsetX,offsetY-Radius,z_left_outer ));

    Set_coordinate_point( 8, D3vector ( offsetX+radius,offsetY,z_right_inner ));
    Set_coordinate_point( 9, D3vector ( offsetX+Radius,offsetY,z_right_outer ));

    Set_coordinate_point( 10, D3vector ( offsetX,offsetY+radius,z_right_inner ));
    Set_coordinate_point( 11, D3vector ( offsetX,offsetY+Radius,z_right_outer ));

    Set_coordinate_point( 12, D3vector ( offsetX-radius,offsetY,z_right_inner ));
    Set_coordinate_point( 13, D3vector ( offsetX-Radius,offsetY,z_right_outer ));

    Set_coordinate_point( 14, D3vector ( offsetX,offsetY-radius,z_right_inner ));
    Set_coordinate_point( 15, D3vector ( offsetX,offsetY-Radius,z_right_outer ));

    int num_blocks = 1;
    for ( int i = 0; i < num_blocks; ++i ) {
      Set_hexahedron ( i, i*num_points_per_face + 0, i*num_points_per_face + 2,
                          i*num_points_per_face + 6, i*num_points_per_face + 4,
                          ( i+1 ) * num_points_per_face + 0, ( i+1 ) * num_points_per_face +2,
                          ( i+1 ) * num_points_per_face + 6, ( i+1 ) * num_points_per_face +4 );
    }
    for ( int j = 0; j < num_points_per_face; j=j+2) {
      for ( int i = 0; i < num_blocks; ++i) {
        int myIndex = (j/2+1)*num_blocks +i;
        Set_hexahedron ( myIndex,
                         i*num_points_per_face     +  j+0,
                         i*num_points_per_face     +  j+1,
                         i*num_points_per_face     +((j+2) % num_points_per_face),
                         i*num_points_per_face     +((j+3) % num_points_per_face),
                         (i+1)*num_points_per_face +  j+0,
                         (i+1)*num_points_per_face +  j+1,
                         (i+1)*num_points_per_face +((j+2) % num_points_per_face),
                         (i+1)*num_points_per_face +((j+3) % num_points_per_face)
                       );
     }
    }
  construction_hexahedron_points_done();


    //inner block: corner ids: 0 2 4 6 8 10 12 14

  Set_transformation_face(0,2,4,6,transform_left_lens_inner_quad); // z=0 plane transform_left_lens_inner_quad
  Set_transformation_face(0,2,8,10,transform_inner_faces_NE_quad);  //transform_inner_quad_NE transform_inner_quad_NW transform_inner_quad_SW transform_inner_quad_SE
  Set_transformation_face(2,4,10,12,transform_inner_faces_NE_quad); //??
  Set_transformation_face(4,6,12,14,transform_inner_faces_NE_quad); //??
  Set_transformation_face(0,6,8,14,transform_inner_faces_NE_quad); //??
  Set_transformation_face(8,10,12,14,transform_right_lens_inner_quad); // z = right plane transform_right_lens_inner_quad

  //top-east block: corner ids: 0 1 2 3 8 9 10 11
  Set_transformation_face(0,1,2,3,transform_left_lens_diag_NE_quad);
  Set_transformation_face(0,1,8,9,transform_diag_inner_faces_NE_quad);
  Set_transformation_face(0,2,8,10,transform_inner_faces_NE_quad);
  Set_transformation_face(2,3,10,11,transform_diag_inner_faces_NE_quad);
  Set_transformation_face(1,3,9,11,transform_outer_boundary_NE); //
  Set_transformation_face(8,9,10,11,transform_right_lens_diag_NE_quad);

  //north-west block: corner ids: 2 3 4 5 10 11 12 13
  Set_transformation_face(2,3,4,5,transform_left_lens_diag_NW_quad);//transform_left_lens_diag_NW_quad
  Set_transformation_face(2,3,10,11,transform_diag_inner_faces_NE_quad);
  Set_transformation_face(2,4,10,12,transform_inner_faces_NE_quad);
  Set_transformation_face(4,5,12,13,transform_diag_inner_faces_NE_quad);
  Set_transformation_face(3,5,11,13,transform_outer_boundary_NW); //transform_outer_boundary_NW
  Set_transformation_face(10,11,12,13,transform_right_lens_diag_NW_quad); //transform_right_lens_diag_NW_quad

  //bottom-west block: corner ids: 4 5 6 7 12 13 14 15
  //eher: S-W-BLOCK
  Set_transformation_face(4,5,6,7,transform_left_lens_diag_SW_quad);
  Set_transformation_face(4,5,12,13,transform_diag_inner_faces_NE_quad);
  Set_transformation_face(4,6,12,14,transform_inner_faces_NE_quad);
  Set_transformation_face(6,7,14,15,transform_diag_inner_faces_NE_quad);
  Set_transformation_face(5,7,13,15,transform_outer_boundary_SW); //
  Set_transformation_face(12,13,14,15,transform_right_lens_diag_SW_quad);

  //bottom-east block: corner ids: 0 1 6 7 8 9 14 15
  //eher: S-E-BLOCK
  Set_transformation_face(0,1,6,7,transform_left_lens_diag_SE_quad);
  Set_transformation_face(0,1,8,9,transform_diag_inner_faces_NE_quad);
  Set_transformation_face(0,6,8,14,transform_inner_faces_NE_quad);
  Set_transformation_face(6,7,14,15,transform_diag_inner_faces_NE_quad);
  Set_transformation_face(1,7,9,15,transform_outer_boundary_SE);
  Set_transformation_face(8,9,14,15,transform_right_lens_diag_SE_quad);



  construction_done();
}

Lens_Geometry_cutted_edges::Lens_Geometry_cutted_edges(double RadiusLeft, double RadiusRight, double MechanicalRadiusLeft,double MechanicalRadiusRight, double thickness, double curvatureLeft, double curvatureRight, double offsetX, double offsetY, double offsetZ, bool inner_grid_arched, double radius)
{
    double maxRadius = std::max(std::max(RadiusLeft,RadiusRight),std::max(MechanicalRadiusLeft,MechanicalRadiusRight));

    bool cutLeft = false;
    if (RadiusLeft < MechanicalRadiusLeft)
    {
        cutLeft = true;
    }
    bool cutRight = false;
    if (RadiusRight < MechanicalRadiusRight)
    {
        cutRight = true;
    }
    if (MechanicalRadiusLeft < RadiusLeft)
    {
        std::cout << "warning : mech radius < radius\n set to radiusLeft";
        MechanicalRadiusLeft = RadiusLeft;
        MechanicalRadiusLeft = maxRadius;
    }
    if (MechanicalRadiusRight < RadiusRight)
    {
        std::cout << "warning : mech radius < radius\n set to radiusRight";
        MechanicalRadiusRight = RadiusRight;
        MechanicalRadiusRight = maxRadius;
    }
    if (RadiusLeft != RadiusRight)
    {
        std::cout << "no support yet for RadiusLeft != RadiusRight ! RadiusLeft will be used for both\n";
    }
    double Radius = RadiusLeft;
    radius = RadiusLeft / 2;


    //assert(Radius > radius);
    size_pointer_global_data = 14;
    if (fabs(curvatureLeft) >1e10)
    {
        curvatureLeft = 1e10;
    }
    if (fabs(curvatureRight) >1e10)
    {
        curvatureRight = 1e10;
    }

    double inversFocal = (1.5 - 1.0) * ( (1.0 / curvatureLeft) - ( 1.0 / -curvatureRight ) + ( 1.5 - 1.0) * thickness / ( 1.5 * -curvatureRight * curvatureLeft ));

    focalLength_ = 1.0 / inversFocal;
    //cout << "Focal length of Lens is " << focalLength_ << " (assume n = 1.5 ) ." << endl;


    thickness_ = thickness;
    offsetX_ = offsetX;
    offsetY_ = offsetY;
    offsetZ_ = offsetZ;

    pointer_global_data = new double[size_pointer_global_data];
    pointer_global_data[0] = Radius;
    pointer_global_data[1] = radius;
    pointer_global_data[2] = curvatureLeft;
    pointer_global_data[3] = curvatureRight;
    pointer_global_data[4] = thickness;
    pointer_global_data[5] = offsetZ;

    assert ((thickness > 0.0) &&  "Thickness must be positive.");
    assert ((Radius > 0.0) &&  "Radius must be positive.");
    assert ((radius > 0.0) &&  "radius must be positive.");



    constructionParameters.resize(12);
    //obacht : correct? should be! not enetiryls sure
    constructionParameters[0] = Radius;
    constructionParameters[0] = maxRadius;
    constructionParameters[1] = radius;
    constructionParameters[2] = thickness;
    constructionParameters[3] = curvatureLeft;
    constructionParameters[4] = curvatureRight;
    constructionParameters[5] = offsetX;
    constructionParameters[6] = offsetY;
    constructionParameters[7] = offsetZ;
    constructionBoolArched = inner_grid_arched;
    // maybe adjust, if lens is negative (most thick part not in the centre but at the edges!)
    relativeCoordVector.push_back(0.0);
    relativeCoordVector.push_back(thickness);
    typFuerSlice = lensUgTyp;

    int num_blocks_z_dir = 1;
    assert(num_blocks_z_dir>0);
    int num_blocks_total = 9 * num_blocks_z_dir;
    int num_points_per_face = 8;
    int numPoints = num_points_per_face*(num_blocks_z_dir+1) ;
    numPoints = 24;

    double z_left_inner  = -sqrt( pow(curvatureLeft ,2) - radius * radius ) + curvatureLeft  + offsetZ;
    double z_left_outer  = -sqrt( pow(curvatureLeft ,2) - Radius * Radius ) + curvatureLeft  + offsetZ;
    double z_right_inner = +sqrt( pow(curvatureRight,2) - radius * radius ) - curvatureRight + offsetZ + thickness;
    double z_right_outer = +sqrt( pow(curvatureRight,2) - Radius * Radius ) - curvatureRight + offsetZ + thickness;

    double z_left_edge;
    if (cutLeft)
    {        z_left_edge = z_left_outer;    }
    else
    {        z_left_edge = -sqrt( pow(curvatureLeft ,2) - MechanicalRadiusLeft * MechanicalRadiusLeft ) + curvatureLeft  + offsetZ;
    }

    double z_right_edge;
    if (cutRight)
    {        z_right_edge = z_right_outer;    }
    else
    {        z_right_edge = +sqrt( pow(curvatureRight,2) - MechanicalRadiusRight * MechanicalRadiusRight ) - curvatureRight + offsetZ + thickness;    }

    //for testing purpose : right edge is calculated, such that Radius defines the z-plane, at which the right lens side is cut.


    if (curvatureLeft< 0.0)
    {
        z_left_inner  = +sqrt( pow(curvatureLeft ,2) - radius * radius ) + curvatureLeft  + offsetZ;
        z_left_outer  = +sqrt( pow(curvatureLeft ,2) - Radius * Radius ) + curvatureLeft  + offsetZ;

        if (cutLeft)
        {        z_left_edge = z_left_outer;    }
        else
        {        z_left_edge = +sqrt( pow(curvatureLeft ,2) - MechanicalRadiusLeft * MechanicalRadiusLeft ) + curvatureLeft  + offsetZ;
        }
      }
    if (curvatureRight< 0.0)
    {
        z_right_inner = -sqrt( pow(curvatureRight,2) - radius * radius ) - curvatureRight + offsetZ + thickness;
        z_right_outer = -sqrt( pow(curvatureRight,2) - Radius * Radius ) - curvatureRight + offsetZ + thickness;
        if (cutRight)
        {        z_right_edge = z_right_outer;    }
        else
        {        z_right_edge = -sqrt( pow(curvatureRight,2) - MechanicalRadiusRight * MechanicalRadiusRight ) - curvatureRight + offsetZ + thickness;    }

    }

    if (curvatureLeft == 0.0)
    {
        z_left_inner  = 0.0  + offsetZ;
        z_left_outer  = 0.0  + offsetZ;
        z_left_edge  = 0.0  + offsetZ;
    }
    if (curvatureRight == 0.0)
    {
        z_right_inner = 0.0 + offsetZ + thickness;
        z_right_outer = 0.0 + offsetZ + thickness;
        z_right_edge  = 0.0 + offsetZ + thickness;
    }


    if (cutLeft)
    {

        constructionParameters[8] = z_left_edge ;
    }
    else
    {
        constructionParameters[8] = std::min(offsetZ,z_left_outer);
    }
    if (cutRight)
    {

        constructionParameters[9] = z_right_edge ;
    }
    else
    {
        constructionParameters[9] = std::max(offsetZ + thickness,z_right_outer);
    }


    if (cutLeft)
    {constructionParameters[10] = Radius;}
    else
    {
       constructionParameters[10] = 0;
    }
    if (cutRight)
    {constructionParameters[11] = Radius;}
    else
    {
       constructionParameters[11] = 0;
    }


//    if (curvatureRight < 0 && curvatureLeft < 0)
//    {
//        assert (((fabs(z_right_outer - z_left_outer) - thickness) < 0.0) &&  "Lens has negative thickness at edge! Make lens more thick or increase radius of curvature.");

//    }
    if (curvatureRight > 0 && curvatureLeft > 0)
    {
        assert ((((z_right_outer - z_left_outer) ) > 0.0) &&  "Lens has negative thickness at edge! Make lens more thick or increase radius of curvature.");

    }


    pointer_global_data[6] = z_left_inner;
    pointer_global_data[7] = z_left_outer;
    pointer_global_data[8] = z_right_inner;
    pointer_global_data[9] = z_right_outer;
    pointer_global_data[10] = z_right_edge;
    pointer_global_data[11] = z_left_edge;
    pointer_global_data[12] = MechanicalRadiusRight;
    pointer_global_data[13] = MechanicalRadiusLeft;
    boundaryPoints.push_back(D3vector{-MechanicalRadiusRight-1,-MechanicalRadiusRight-1,constructionParameters[8]-1});
    boundaryPoints.push_back(D3vector{+MechanicalRadiusRight+1,+MechanicalRadiusRight+1,constructionParameters[9]+1});

//    cout << "fixed lens param., channge again " << endl;
//    radius = 1;
//    Radius = 2;
//    z_left_inner = z_left_outer = 0;
//    z_right_inner = z_right_outer = 1;



    Set_number_points ( numPoints );
    Set_number_hexahedra ( num_blocks_total );


    Set_coordinate_point( 0, D3vector ( offsetX+radius,offsetY,z_left_inner ));
    Set_coordinate_point( 1, D3vector ( offsetX+Radius,offsetY,z_left_outer ));

    Set_coordinate_point( 2, D3vector ( offsetX       ,offsetY+radius,z_left_inner ));
    Set_coordinate_point( 3, D3vector ( offsetX       ,offsetY+Radius,z_left_outer ));

    Set_coordinate_point( 4, D3vector ( offsetX-radius,offsetY,z_left_inner ));
    Set_coordinate_point( 5, D3vector ( offsetX-Radius,offsetY,z_left_outer ));

    Set_coordinate_point( 6, D3vector ( offsetX       ,offsetY-radius,z_left_inner ));
    Set_coordinate_point( 7, D3vector ( offsetX       ,offsetY-Radius,z_left_outer ));




    Set_coordinate_point( 8, D3vector ( offsetX+radius,offsetY,z_right_inner ));
    Set_coordinate_point( 9, D3vector ( offsetX+Radius,offsetY,z_right_outer ));

    Set_coordinate_point( 10, D3vector ( offsetX        ,offsetY+radius,z_right_inner ));
    Set_coordinate_point( 11, D3vector ( offsetX        ,offsetY+Radius,z_right_outer ));

    Set_coordinate_point( 12, D3vector ( offsetX-radius,offsetY,z_right_inner ));
    Set_coordinate_point( 13, D3vector ( offsetX-Radius,offsetY,z_right_outer ));

    Set_coordinate_point( 14, D3vector ( offsetX        ,offsetY-radius,z_right_inner ));
    Set_coordinate_point( 15, D3vector ( offsetX        ,offsetY-Radius,z_right_outer ));


    Set_coordinate_point( 16, D3vector ( offsetX+MechanicalRadiusLeft,offsetY,z_left_edge ));
    Set_coordinate_point( 17, D3vector ( offsetX,offsetY+MechanicalRadiusLeft,z_left_edge ));

    Set_coordinate_point( 18, D3vector ( offsetX-MechanicalRadiusLeft,offsetY,z_left_edge ));
    Set_coordinate_point( 19, D3vector ( offsetX,offsetY-MechanicalRadiusLeft,z_left_edge ));

    Set_coordinate_point( 20, D3vector ( offsetX+MechanicalRadiusRight,offsetY,z_right_edge ));
    Set_coordinate_point( 21, D3vector ( offsetX,offsetY+MechanicalRadiusRight,z_right_edge ));

    Set_coordinate_point( 22, D3vector ( offsetX-MechanicalRadiusRight,offsetY,z_right_edge ));
    Set_coordinate_point( 23, D3vector ( offsetX,offsetY-MechanicalRadiusRight,z_right_edge ));



    int num_blocks = 1;
    for ( int i = 0; i < num_blocks; ++i ) {
      Set_hexahedron ( i, i*num_points_per_face + 0, i*num_points_per_face + 2,
                          i*num_points_per_face + 6, i*num_points_per_face + 4,
                          ( i+1 ) * num_points_per_face + 0, ( i+1 ) * num_points_per_face +2,
                          ( i+1 ) * num_points_per_face + 6, ( i+1 ) * num_points_per_face +4 );
    }
    //outer hexhedrons
    for ( int j = 0; j < num_points_per_face; j=j+2) {
      for ( int i = 0; i < num_blocks; ++i) {
        int myIndex = (j/2+1)*num_blocks +i;
        Set_hexahedron ( myIndex,
                         i*num_points_per_face     +  j+0,
                         i*num_points_per_face     +  j+1,
                         i*num_points_per_face     +((j+2) % num_points_per_face),
                         i*num_points_per_face     +((j+3) % num_points_per_face),
                         (i+1)*num_points_per_face +  j+0,
                         (i+1)*num_points_per_face +  j+1,
                         (i+1)*num_points_per_face +((j+2) % num_points_per_face),
                         (i+1)*num_points_per_face +((j+3) % num_points_per_face)
                       );
     }
    }
    std::cout << "outer blocks " << std::endl;
    //outer hexhedrons

        Set_hexahedron ( 5, 1,16,3,17,1+8,16+4,3+8,17+4);
        Set_hexahedron ( 6, 3,17,5,18,3+8,17+4,5+8,18+4);
        Set_hexahedron ( 7, 5,18,7,19,5+8,18+4,7+8,19+4);
        Set_hexahedron ( 8, 7,19,1,16,7+8,19+4,1+8,16+4);


  construction_hexahedron_points_done();

//put to back again

  //top-east outer block: corner ids: 0 1  2  3 8 9  10 11
  //top-east outer block: corner ids: 1 16 3 17 9 20 11 21
  //1,16,3,17,1+8,16+4,3+8,17+4
  Set_transformation_face(1,16,3,17,transform_left_lens_diag_NE_quad_cut);
  Set_transformation_face(1,16,9,20,transform_diag_inner_faces_NE_quad_cut);
  Set_transformation_face(1,3,9,11,transform_outer_boundary_NE); // inner face :
  Set_transformation_face(3,17,11,21,transform_diag_inner_faces_NE_quad_cut);
  Set_transformation_face(16,17,20,21,transform_outer_boundary_NE_cut); //
  Set_transformation_face(9,20,11,21,transform_right_lens_diag_NE_quad_cut);

  //north-west outer block: corner ids: 2 3  4 5  10 11 12 13
  //north-west outer block: corner ids: 3 17 5 18 11 21 13 22
  //3,17,5,18,3+8,17+4,5+8,18+4
  Set_transformation_face(3,17,5,18,transform_left_lens_diag_NW_quad_cut);//transform_left_lens_diag_NW_quad
  Set_transformation_face(3,17,11,21,transform_diag_inner_faces_NE_quad_cut);
  Set_transformation_face(3,5,11,13,transform_outer_boundary_NW); // inner face : no transform necessary
  Set_transformation_face(5,18,13,22,transform_diag_inner_faces_NE_quad_cut);
  Set_transformation_face(17,18,21,22,transform_outer_boundary_NW_cut); //transform_outer_boundary_NW
  Set_transformation_face(11,21,13,22,transform_right_lens_diag_NW_quad_cut); //transform_right_lens_diag_NW_quad

  //bottom-west outer block: corner ids: 4 5  6 7  12 13 14 15
  //bottom-west outer block: corner ids: 5 18 7 19 13 22 15 23
  //5,18,7,19,5+8,18+4,7+8,19+4
  Set_transformation_face(5,18,7,19,transform_left_lens_diag_SW_quad_cut);
  Set_transformation_face(5,18,13,22,transform_diag_inner_faces_NE_quad_cut);
  Set_transformation_face(5,7,13,15,transform_outer_boundary_SW); // inner face : no transform necessary
  Set_transformation_face(7,19,15,23,transform_diag_inner_faces_NE_quad_cut);
  Set_transformation_face(18,19,22,23,transform_outer_boundary_SW_cut); //
  Set_transformation_face(13,22,15,23,transform_right_lens_diag_SW_quad_cut);

  //bottom-east outer block: corner ids: 0 1  6 7  8  9 14 15
  //bottom-east outer block: corner ids: 1 16 7 19 9 20 15 23
  //7,19,1,16,7+8,19+4,1+8,16+4
  Set_transformation_face(1,16,7,19,transform_left_lens_diag_SE_quad_cut);
  Set_transformation_face(1,16,9,20,transform_diag_inner_faces_NE_quad_cut);
  Set_transformation_face(1,7,9,15,transform_outer_boundary_SE); // inner face : no transform necessary
  Set_transformation_face(7,19,15,23,transform_diag_inner_faces_NE_quad_cut);
  Set_transformation_face(16,19,20,23,transform_outer_boundary_SE_cut);
  Set_transformation_face(9,20,15,23,transform_right_lens_diag_SE_quad_cut);

    //inner block: corner ids: 0 2 4 6 8 10 12 14

  Set_transformation_face(0,2,4,6,transform_left_lens_inner_quad); // z=0 plane transform_left_lens_inner_quad
  Set_transformation_face(0,2,8,10,transform_inner_faces_NE_quad);  //transform_inner_quad_NE transform_inner_quad_NW transform_inner_quad_SW transform_inner_quad_SE
  Set_transformation_face(2,4,10,12,transform_inner_faces_NE_quad); //??
  Set_transformation_face(4,6,12,14,transform_inner_faces_NE_quad); //??
  Set_transformation_face(0,6,8,14,transform_inner_faces_NE_quad); //??
  Set_transformation_face(8,10,12,14,transform_right_lens_inner_quad); // z = right plane transform_right_lens_inner_quad

  //top-east block: corner ids: 0 1 2 3 8 9 10 11
  Set_transformation_face(0,1,2,3,transform_left_lens_diag_NE_quad);
  Set_transformation_face(0,1,8,9,transform_diag_inner_faces_NE_quad);
  Set_transformation_face(0,2,8,10,transform_inner_faces_NE_quad);
  Set_transformation_face(2,3,10,11,transform_diag_inner_faces_NE_quad);
  Set_transformation_face(1,3,9,11,transform_outer_boundary_NE); //
  Set_transformation_face(8,9,10,11,transform_right_lens_diag_NE_quad);


  //north-west block: corner ids: 2 3 4 5 10 11 12 13
  Set_transformation_face(2,3,4,5,transform_left_lens_diag_NW_quad);//transform_left_lens_diag_NW_quad
  Set_transformation_face(2,3,10,11,transform_diag_inner_faces_NE_quad);
  Set_transformation_face(2,4,10,12,transform_inner_faces_NE_quad);
  Set_transformation_face(4,5,12,13,transform_diag_inner_faces_NE_quad);
  Set_transformation_face(3,5,11,13,transform_outer_boundary_NW); //transform_outer_boundary_NW
  Set_transformation_face(10,11,12,13,transform_right_lens_diag_NW_quad); //transform_right_lens_diag_NW_quad

  //bottom-west block: corner ids: 4 5 6 7 12 13 14 15
  //eher: S-W-BLOCK
  Set_transformation_face(4,5,6,7,transform_left_lens_diag_SW_quad);
  Set_transformation_face(4,5,12,13,transform_diag_inner_faces_NE_quad);
  Set_transformation_face(4,6,12,14,transform_inner_faces_NE_quad);
  Set_transformation_face(6,7,14,15,transform_diag_inner_faces_NE_quad);
  Set_transformation_face(5,7,13,15,transform_outer_boundary_SW); //
  Set_transformation_face(12,13,14,15,transform_right_lens_diag_SW_quad);

  //bottom-east block: corner ids: 0 1 6 7 8 9 14 15
  //eher: S-E-BLOCK
  Set_transformation_face(0,1,6,7,transform_left_lens_diag_SE_quad);
  Set_transformation_face(0,1,8,9,transform_diag_inner_faces_NE_quad);
  Set_transformation_face(0,6,8,14,transform_inner_faces_NE_quad);
  Set_transformation_face(6,7,14,15,transform_diag_inner_faces_NE_quad);
  Set_transformation_face(1,7,9,15,transform_outer_boundary_SE);
  Set_transformation_face(8,9,14,15,transform_right_lens_diag_SE_quad);




  // WSDdir3D, ESDdir3D, WNDdir3D, ENDdir3D, WSTdir3D, ESTdir3D, WNTdir3D, ENTdir3D

  construction_done();

}