examples_ug_optics.cc

#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();
}