Vector Quadrupole::B(const Vector &position) {
    
    double dX,dY,dZ; //offsets from quadrupole center
    double Bx_,By_,Bz_; //magnetic field in the transformed coordinates
    
    //Calculate offset
    dX=position.x-center.x;
    dY=position.y-center.y;
    dZ=position.z-center.z;

    //Transform 
    double Px=dX*cos(Beta)+dZ*sin(Beta);
    double Py=-dX*sin(Alpha)*sin(Beta)+dY*cos(Alpha)+dZ*sin(Alpha)*cos(Beta);
    double Pz=-dX*cos(Alpha)*sin(Beta)-dY*sin(Alpha)+dZ*cos(Alpha)*cos(Beta);

    dX=Px*cos(rot)+Py*sin(rot);
    dY=-Px*sin(rot)+Py*cos(rot);
    Px=dX;
    Py=dY;

    if (Pz>=0.0 && Pz<=L) { //if inside the quadrupole
        Bx_=-K*Py;
        By_=-K*Px;
    } else {
        Bx_=0.0;
        By_=0.0;
    }
    Bz_=0.0;

    //Inverse transformation
    Px=Bx_*cos(rot)-By_*sin(rot);
    Py=Bx_*sin(rot)+By_*cos(rot);
    Pz=Bz_;

    return Vector(Px*cos(Beta)-Py*sin(Alpha)*sin(Beta)-Pz*cos(Alpha)*sin(Beta),
        Py*cos(Alpha) -Pz*sin(Alpha),
        Px*sin(Beta)+Py*sin(Alpha)*cos(Beta) +Pz*cos(Alpha)*cos(Beta);
}