Skip to content

Commit

Permalink
curveplanarfourier.cpp: Reduced number of calls to trig functions
Browse files Browse the repository at this point in the history
  • Loading branch information
@landreman
landreman committed Nov 27, 2023
1 parent 73637ae commit 60cb677
Showing 1 changed file with 64 additions and 93 deletions.
157 changes: 64 additions & 93 deletions src/simsoptpp/curveplanarfourier.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,17 +24,13 @@ void CurvePlanarFourier<Array>::gamma_impl(Array& data, Array& quadpoints) {
double phi = 2 * M_PI * quadpoints[k];
double cosphi = cos(phi);
double sinphi = sin(phi);
double cosiphi = 0;
double siniphi = 0;
for (int i = 0; i < order+1; ++i) {
cosiphi = cos(i*phi);
data(k, 0) += rc[i] * cosiphi * cosphi;
data(k, 1) += rc[i] * cosiphi * sinphi;
}
data(k, 0) = rc[0] * cosphi;
data(k, 1) = rc[0] * sinphi;
for (int i = 1; i < order+1; ++i) {
cosiphi = cos(i*phi);
siniphi = sin(i*phi);
data(k, 0) += rs[i-1] * siniphi * cosphi;
data(k, 1) += rs[i-1] * siniphi * sinphi;
data(k, 0) += (rc[i] * cosiphi + rs[i-1] * siniphi) * cosphi;
data(k, 1) += (rc[i] * cosiphi + rs[i-1] * siniphi) * sinphi;
}
}
for (int m = 0; m < numquadpoints; ++m) {
Expand All @@ -56,23 +52,20 @@ void CurvePlanarFourier<Array>::gammadash_impl(Array& data) {
double inv_sqrt_s = inv_magnitude();
Array q_norm = q * inv_sqrt_s;

double cosiphi, siniphi;
for (int k = 0; k < numquadpoints; ++k) {
double phi = 2 * M_PI * quadpoints[k];
double cosphi = cos(phi);
double sinphi = sin(phi);
double cosiphi = 0;
double siniphi = 0;
for (int i = 0; i < order+1; ++i) {
cosiphi = cos(i*phi);
siniphi = sin(i*phi);
data(k, 0) += rc[i] * ( -(i) * siniphi * cosphi - cosiphi * sinphi);
data(k, 1) += rc[i] * ( -(i) * siniphi * sinphi + cosiphi * cosphi);
}
data(k, 0) = rc[0] * (-sinphi);
data(k, 1) = rc[0] * (cosphi);
for (int i = 1; i < order+1; ++i) {
cosiphi = cos(i*phi);
siniphi = sin(i*phi);
data(k, 0) += rs[i-1] * ( (i) * cosiphi * cosphi - siniphi * sinphi);
data(k, 1) += rs[i-1] * ( (i) * cosiphi * sinphi + siniphi * cosphi);
data(k, 0) += rc[i] * ( -(i) * siniphi * cosphi - cosiphi * sinphi)
+ rs[i-1] * ( (i) * cosiphi * cosphi - siniphi * sinphi);
data(k, 1) += rc[i] * ( -(i) * siniphi * sinphi + cosiphi * cosphi)
+ rs[i-1] * ( (i) * cosiphi * sinphi + siniphi * cosphi);
}
}

Expand All @@ -96,23 +89,20 @@ void CurvePlanarFourier<Array>::gammadashdash_impl(Array& data) {

Array q_norm = q * inv_magnitude();

double cosiphi, siniphi;
for (int k = 0; k < numquadpoints; ++k) {
double phi = 2 * M_PI * quadpoints[k];
double cosphi = cos(phi);
double sinphi = sin(phi);
double cosiphi = 0;
double siniphi = 0;
for (int i = 0; i < order+1; ++i) {
cosiphi = cos(i*phi);
siniphi = sin(i*phi);
data(k, 0) += rc[i] * (+2*(i)*siniphi*sinphi-(i*i+1)*cosiphi*cosphi);
data(k, 1) += rc[i] * (-2*(i)*siniphi*cosphi-(i*i+1)*cosiphi*sinphi);
}
data(k, 0) = rc[0] * (-cosphi);
data(k, 1) = rc[0] * (-sinphi);
for (int i = 1; i < order+1; ++i) {
cosiphi = cos(i*phi);
siniphi = sin(i*phi);
data(k, 0) += rs[i-1] * (-(i*i+1)*siniphi*cosphi - 2*(i)*cosiphi*sinphi);
data(k, 1) += rs[i-1] * (-(i*i+1)*siniphi*sinphi + 2*(i)*cosiphi*cosphi);
data(k, 0) += rc[i] * (+2*(i)*siniphi*sinphi-(i*i+1)*cosiphi*cosphi)
+ rs[i-1] * (-(i*i+1)*siniphi*cosphi - 2*(i)*cosiphi*sinphi);
data(k, 1) += rc[i] * (-2*(i)*siniphi*cosphi-(i*i+1)*cosiphi*sinphi)
+ rs[i-1] * (-(i*i+1)*siniphi*sinphi + 2*(i)*cosiphi*cosphi);
}
}
data *= 2*M_PI*2*M_PI;
Expand All @@ -135,33 +125,28 @@ void CurvePlanarFourier<Array>::gammadashdashdash_impl(Array& data) {

Array q_norm = q * inv_magnitude();

double cosiphi, siniphi;
for (int k = 0; k < numquadpoints; ++k) {
double phi = 2 * M_PI * quadpoints[k];
double cosphi = cos(phi);
double sinphi = sin(phi);
double cosiphi = 0;
double siniphi = 0;

for (int i = 0; i < order+1; ++i) {
data(k, 0) = rc[0]*(+sinphi);
data(k, 1) = rc[0]*(-cosphi);
for (int i = 1; i < order+1; ++i) {
cosiphi = cos(i*phi);
siniphi = sin(i*phi);
data(k, 0) += rc[i]*(
+(3*i*i + 1)*cosiphi*sinphi
+(i*i + 3)*(i)*siniphi*cosphi
) + rs[i-1]*(
-(i*i+3) * (i) * cosiphi*cosphi
+(3*i*i+1) * siniphi*sinphi
);
data(k, 1) += rc[i]*(
+(i*i + 3)*(i)*siniphi*sinphi
-(3*i*i + 1)*cosiphi*cosphi
);
}
for (int i = 1; i < order+1; ++i) {
cosiphi = cos(i*phi);
siniphi = sin(i*phi);
data(k, 0) += rs[i-1]*(
-(i*i+3) * (i) * cosiphi*cosphi
+(3*i*i+1) * siniphi*sinphi
);
data(k, 1) += rs[i-1]*(
) + rs[i-1]*(
-(i*i+3)*(i)*cosiphi*sinphi
-(3*i*i+1)*siniphi*cosphi
);
Expand All @@ -187,6 +172,7 @@ void CurvePlanarFourier<Array>::dgamma_by_dcoeff_impl(Array& data) {

Array q_norm = q * inv_magnitude();

double cosnphi, sinnphi;
for (int m = 0; m < numquadpoints; ++m) {
double phi = 2 * M_PI * quadpoints[m];
int counter = 0;
Expand All @@ -200,8 +186,9 @@ void CurvePlanarFourier<Array>::dgamma_by_dcoeff_impl(Array& data) {
double sinnphi = 0;

for (int n = 0; n < order+1; ++n) {
i = cos(n*phi) * cosphi;
j = cos(n*phi) * sinphi;
cosnphi = cos(n * phi);
i = cosnphi * cosphi;
j = cosnphi * sinphi;
k = 0;

data(m, 0, counter) = (i - 2 * (q_norm[2] * q_norm[2] + q_norm[3] * q_norm[3]) * i + 2 * (q_norm[1] * q_norm[2] - q_norm[3] * q_norm[0]) * j + 2 * (q_norm[1] * q_norm[3] + q_norm[2] * q_norm[0]) * k);
Expand All @@ -212,8 +199,9 @@ void CurvePlanarFourier<Array>::dgamma_by_dcoeff_impl(Array& data) {
}

for (int n = 1; n < order+1; ++n) {
i = sin(n*phi) * cosphi;
j = sin(n*phi) * sinphi;
sinnphi = sin(n * phi);
i = sinnphi * cosphi;
j = sinnphi * sinphi;
k = 0;

data(m, 0, counter) = (i - 2 * (q_norm[2] * q_norm[2] + q_norm[3] * q_norm[3]) * i + 2 * (q_norm[1] * q_norm[2] - q_norm[3] * q_norm[0]) * j + 2 * (q_norm[1] * q_norm[3] + q_norm[2] * q_norm[0]) * k);
Expand All @@ -223,17 +211,15 @@ void CurvePlanarFourier<Array>::dgamma_by_dcoeff_impl(Array& data) {
counter++;
}

i = 0;
j = 0;
i = rc[0] * cosphi;
j = rc[0] * sinphi;
k = 0;

for (int n = 0; n < order+1; ++n) {
i += rc[n] * cos(n*phi) * cosphi;
j += rc[n] * cos(n*phi) * sinphi;
}
for (int n = 1; n < order+1; ++n) {
i += rs[n-1] * sin(n*phi) * cosphi;
j += rs[n-1] * sin(n*phi) * sinphi;
cosnphi = cos(n * phi);
sinnphi = sin(n * phi);
i += (rc[n] * cosnphi + rs[n-1] * sinnphi) * cosphi;
j += (rc[n] * cosnphi + rs[n-1] * sinnphi) * sinphi;
}
double inv_sqrt_s = inv_magnitude();

Expand Down Expand Up @@ -530,12 +516,11 @@ void CurvePlanarFourier<Array>::dgammadash_by_dcoeff_impl(Array& data) {

Array q_norm = q * inv_magnitude();

double cosnphi, sinnphi;
for (int m = 0; m < numquadpoints; ++m) {
double phi = 2 * M_PI * quadpoints[m];
double cosphi = cos(phi);
double sinphi = sin(phi);
double cosnphi = 0;
double sinnphi = 0;
int counter = 0;
double i;
double j;
Expand Down Expand Up @@ -576,21 +561,19 @@ void CurvePlanarFourier<Array>::dgammadash_by_dcoeff_impl(Array& data) {

}

i = 0;
j = 0;
i = rc[0] * (-sinphi);
j = rc[0] * (cosphi);
k = 0;
for (int n = 0; n < order+1; ++n) {
cosnphi = cos(n*phi);
sinnphi = sin(n*phi);
i += rc[n] * ( -(n) * sinnphi * cosphi - cosnphi * sinphi) * 2 * M_PI;
j += rc[n] * ( -(n) * sinnphi * sinphi + cosnphi * cosphi) * 2 * M_PI;
}
for (int n = 1; n < order+1; ++n) {
cosnphi = cos(n*phi);
sinnphi = sin(n*phi);
i += rs[n-1] * ( (n) * cosnphi * cosphi - sinnphi * sinphi) * 2 * M_PI;
j += rs[n-1] * ( (n) * cosnphi * sinphi + sinnphi * cosphi) * 2 * M_PI;
i += rc[n] * ( -(n) * sinnphi * cosphi - cosnphi * sinphi)
+ rs[n-1] * ( (n) * cosnphi * cosphi - sinnphi * sinphi);
j += rc[n] * ( -(n) * sinnphi * sinphi + cosnphi * cosphi)
+ rs[n-1] * ( (n) * cosnphi * sinphi + sinnphi * cosphi);
}
i *= (2*M_PI);
j *= (2*M_PI);
double inv_sqrt_s = inv_magnitude();

data(m, 0, counter) = (4 * i * (q_norm[2] * q_norm[2] + q_norm[3] * q_norm[3]) * q_norm[0]
Expand Down Expand Up @@ -886,12 +869,11 @@ void CurvePlanarFourier<Array>::dgammadashdash_by_dcoeff_impl(Array& data) {

Array q_norm = q * inv_magnitude();

double cosnphi, sinnphi;
for (int m = 0; m < numquadpoints; ++m) {
double phi = 2 * M_PI * quadpoints[m];
double cosphi = cos(phi);
double sinphi = sin(phi);
double cosnphi = 0;
double sinnphi = 0;
int counter = 0;
double i;
double j;
Expand Down Expand Up @@ -932,20 +914,16 @@ void CurvePlanarFourier<Array>::dgammadashdash_by_dcoeff_impl(Array& data) {
counter++;
}

i = 0;
j = 0;
i = rc[0] * (-cosphi);
j = rc[0] * (-sinphi);
k = 0;
for (int n = 0; n < order+1; ++n) {
cosnphi = cos(n*phi);
sinnphi = sin(n*phi);
i += rc[n] * (+2*(n)*sinnphi*sinphi-(n*n+1)*cosnphi*cosphi);
j += rc[n] * (-2*(n)*sinnphi*cosphi-(n*n+1)*cosnphi*sinphi);
}
for (int n = 1; n < order+1; ++n) {
cosnphi = cos(n*phi);
sinnphi = sin(n*phi);
i += rs[n-1] * (-(n*n+1)*sinnphi*cosphi - 2*(n)*cosnphi*sinphi);
j += rs[n-1] * (-(n*n+1)*sinnphi*sinphi + 2*(n)*cosnphi*cosphi);
i += rc[n] * (+2*(n)*sinnphi*sinphi-(n*n+1)*cosnphi*cosphi)
+ rs[n-1] * (-(n*n+1)*sinnphi*cosphi - 2*(n)*cosnphi*sinphi);
j += rc[n] * (-2*(n)*sinnphi*cosphi-(n*n+1)*cosnphi*sinphi)
+ rs[n-1] * (-(n*n+1)*sinnphi*sinphi + 2*(n)*cosnphi*cosphi);
}
i *= 2*M_PI*2*M_PI;
j *= 2*M_PI*2*M_PI;
Expand Down Expand Up @@ -1244,12 +1222,11 @@ void CurvePlanarFourier<Array>::dgammadashdashdash_by_dcoeff_impl(Array& data) {

Array q_norm = q * inv_magnitude();

double cosnphi, sinnphi;
for (int m = 0; m < numquadpoints; ++m) {
double phi = 2 * M_PI * quadpoints[m];
double cosphi = cos(phi);
double sinphi = sin(phi);
double cosnphi = 0;
double sinnphi = 0;
int counter = 0;
double i;
double j;
Expand Down Expand Up @@ -1302,29 +1279,23 @@ void CurvePlanarFourier<Array>::dgammadashdashdash_by_dcoeff_impl(Array& data) {
counter++;
}

i = 0;
j = 0;
i = rc[0]*(sinphi);
j = rc[0]*(-cosphi);
k = 0;
for (int n = 0; n < order+1; ++n) {
for (int n = 1; n < order+1; ++n) {
cosnphi = cos(n*phi);
sinnphi = sin(n*phi);
i += rc[n]*(
+(3*n*n + 1)*cosnphi*sinphi
+(n*n + 3)*(n)*sinnphi*cosphi
) + rs[n-1]*(
-(n*n+3) * (n) * cosnphi*cosphi
+(3*n*n+1) * sinnphi*sinphi
);
j += rc[n]*(
+(n*n + 3)*(n)*sinnphi*sinphi
-(3*n*n + 1)*cosnphi*cosphi
);
}
for (int n = 1; n < order+1; ++n) {
cosnphi = cos(n*phi);
sinnphi = sin(n*phi);
i += rs[n-1]*(
-(n*n+3) * (n) * cosnphi*cosphi
+(3*n*n+1) * sinnphi*sinphi
);
j += rs[n-1]*(
) + rs[n-1]*(
-(n*n+3)*(n)*cosnphi*sinphi
-(3*n*n+1)*sinnphi*cosphi
);
Expand Down

0 comments on commit 60cb677

Please sign in to comment.