``````#ifndef RTS_QUAD_H

//enable CUDA_CALLABLE macro
#include <stim/cuda/callable.h>
#include <stim/math/vector.h>
#include <stim/math/triangle.h>
#include <stim/math/quaternion.h>
#include <iostream>
#include <iomanip>
#include <algorithm>

namespace stim{

//template for a quadangle class in ND space
template <class T, int N = 3>
{
/*
B------------------>C
^                   ^
|                   |
Y                   |
|                   |
|                   |
A---------X-------->O
*/

/*T A[N];
T B[N];
T C[N];*/

rts::vec<T, N> A;
rts::vec<T, N> X;
rts::vec<T, N> Y;

{

}

CUDA_CALLABLE quad(vec<T, N> a, vec<T, N> b, vec<T, N> c)
{

A = a;
Y = b - a;
X = c - a - Y;

}

/*******************************************************************
Constructor - create a quad from a position, normal, and rotation
*******************************************************************/
CUDA_CALLABLE quad(rts::vec<T, N> c, rts::vec<T, N> normal, T width, T height, T theta)
{

//compute the X direction - start along world-space X
Y = rts::vec<T, N>(0, 1, 0);
if(Y == normal)
Y = rts::vec<T, N>(0, 0, 1);

X = Y.cross(normal).norm();

std::cout<<X<<std::endl;

//rotate the X axis by theta radians
rts::quaternion<T> q;
q.CreateRotation(theta, normal);
X = q.toMatrix3() * X;
Y = normal.cross(X);

//normalize everything
X = X.norm();
Y = Y.norm();

//scale to match the quad width and height
X = X * width;
Y = Y * height;

//set the corner of the plane
A = c - X * 0.5f - Y * 0.5f;

std::cout<<X<<std::endl;
}

//boolean comparison
bool operator==(const quad<T, N> & rhs)
{
if(A == rhs.A && X == rhs.X && Y == rhs.Y)
return true;
else
return false;
}

/*******************************************
Return the normal for the quad
*******************************************/
CUDA_CALLABLE rts::vec<T, N> n()
{
return (X.cross(Y)).norm();
}

CUDA_CALLABLE rts::vec<T, N> p(T a, T b)
{
rts::vec<T, N> result;
//given the two parameters a, b = [0 1], returns the position in world space
result = A + X * a + Y * b;

return result;
}

CUDA_CALLABLE rts::vec<T, N> operator()(T a, T b)
{
return p(a, b);
}

std::string str()
{
std::stringstream ss;

ss<<std::left<<"B="<<setfill('-')<<setw(20)<<A + Y<<">"<<"C="<<A + Y + X<<std::endl;
ss<<setfill(' ')<<setw(23)<<"|"<<"|"<<std::endl<<setw(23)<<"|"<<"|"<<std::endl;
ss<<std::left<<"A="<<setfill('-')<<setw(20)<<A<<">"<<"D="<<A + X;

return ss.str();

}

{
//scales the plane by a scalar value

//compute the center point
rts::vec<T, N> c = A + X*0.5f + Y*0.5f;

result.X = X * rhs;
result.Y = Y * rhs;
result.A = c - result.X*0.5f - result.Y*0.5f;

return result;

}

CUDA_CALLABLE T dist(vec<T, N> p)
{
//compute the distance between a point and this quad

//first break the quad up into two triangles
triangle<T, N> T0(A, A+X, A+Y);
triangle<T, N> T1(A+X+Y, A+X, A+Y);

T d0 = T0.dist(p);
T d1 = T1.dist(p);

if(d0 < d1)
return d0;
else
return d1;
}

CUDA_CALLABLE T dist_max(vec<T, N> p)
{
T da = (A - p).len();
T db = (A+X - p).len();
T dc = (A+Y - p).len();
T dd = (A+X+Y - p).len();

return std::max( da, std::max(db, std::max(dc, dd) ) );
}
};

}	//end namespace rts

template <typename T, int N>
std::ostream& operator<<(std::ostream& os, rts::quad<T, N> R)
{
os<<R.str();
return os;
}

#endif
``````