DESC icon indicating copy to clipboard operation
DESC copied to clipboard

2nd Order NAE Constraints

Open dpanici opened this issue 9 months ago • 2 comments

  • [ ] make tests for this
  • [ ] write out the coefficients in a simpler form (see below)

Simpler way to write coefficients

from qsc import Qsc
import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import fft

name = 'precise QA'
stel = Qsc.from_paper(name)
R0 = stel.R0
Z0 = stel.Z0
dR0 = stel.R0p
ddR0 = stel.R0pp
dZ0 = stel.Z0p
ddZ0 = stel.Z0pp

x1_cos = stel.X1c*stel.normal_cylindrical.transpose()+stel.Y1c*stel.binormal_cylindrical.transpose()
x1_sin = stel.Y1s*stel.binormal_cylindrical.transpose()

x2_0 = stel.X20*stel.normal_cylindrical.transpose()+stel.Y20*stel.binormal_cylindrical.transpose()+stel.Z20*stel.tangent_cylindrical.transpose()
x2_cos = stel.X2c*stel.normal_cylindrical.transpose()+stel.Y2c*stel.binormal_cylindrical.transpose()+stel.Z2c*stel.tangent_cylindrical.transpose()
x2_sin = stel.X2s*stel.normal_cylindrical.transpose()+stel.Y2s*stel.binormal_cylindrical.transpose()+stel.Z2s*stel.tangent_cylindrical.transpose()

X1Rc = x1_cos[0,:]
dX1Rc = np.matmul(stel.d_d_phi, X1Rc)
X1Rs = x1_sin[0,:]
dX1Rs = np.matmul(stel.d_d_phi, X1Rs)
X1thc = x1_cos[1,:]
dX1thc = np.matmul(stel.d_d_phi, X1thc)
X1ths = x1_sin[1,:]
dX1ths = np.matmul(stel.d_d_phi, X1ths)
X1zc = x1_cos[2,:]
dX1zc = np.matmul(stel.d_d_phi, X1zc)
X1zs = x1_sin[2,:]
dX1zs = np.matmul(stel.d_d_phi, X1zs)
X2R0 = x2_0[0,:]
X2th0 = x2_0[1,:]
X2z0 = x2_0[2,:]
X2Rs = x2_sin[0,:]
X2ths = x2_sin[1,:]
X2zs = x2_sin[2,:]
X2Rc = x2_cos[0,:]
X2thc = x2_cos[1,:]
X2zc = x2_cos[2,:]

R1c = X1Rc - dR0/R0*X1thc
R1s = X1Rs - dR0/R0*X1ths

Z1c = X1zc - dZ0/R0*X1thc
Z1s = X1zs - dZ0/R0*X1ths

R20 = X2R0 - dR0**2*(X1thc**2+X1ths**2)/2/R0**3 + 1/4/R0*(X1thc**2 + X1ths**2 - 4*X2th0*dR0 - 2*X1thc*dX1Rc - 2*X1ths*dX1Rs) + \
    1/4/R0**2*(2*X1Rc*X1thc*dR0 + 2*X1Rs*X1ths*dR0 + 2*X1thc*dR0*dX1thc + 2*X1ths*dR0*dX1ths + 3*(X1thc**2+X1ths**2)*ddR0)
R2s = X2Rs - X1thc*X1ths*dR0**2/R0**3 -(2*X2ths*dR0 + X1ths*dX1Rc + X1thc*(-X1ths + dX1Rs))/(2*R0) + \
    (X1Rs*X1thc*dR0 + X1Rc*X1ths*dR0 + X1ths*dR0*dX1thc + X1thc*dR0*dX1ths + X1thc*X1ths*ddR0)/(2*R0**2)
R2c = X2Rc + (-X1thc**2 + X1ths**2)*dR0**2/(2*R0**3)+1/4/R0*(X1thc**2 - X1ths**2 - 4*X2thc*dR0 - \
    2*X1thc*dX1Rc + 2*X1ths*dX1Rs) + 1/4/R0**2*(2*X1Rc*X1thc*dR0 - 2*X1Rs*X1ths*dR0 + 2*X1thc*dR0*dX1thc - \
    2*X1ths*dR0*dX1ths + (X1thc**2-X1ths**2)*ddR0)

Z20 = X2z0 - (X1thc*dX1zc + X1ths*dX1zs + 2*X2th0*dZ0)/2/R0 + \
    (X1Rc*X1thc + X1Rs*X1ths + X1thc*dX1thc + X1ths*dX1ths)*dZ0/2/R0**2 + \
    (X1thc**2 + X1ths**2)*ddZ0/4/R0**2 - (X1thc**2 + X1ths**2)*dR0*dZ0/2/R0**3 

Z2s = X2zs - (X1ths*dX1zc + X1thc*dX1zs + 2*X2ths*dZ0)/2/R0 + \
    ((X1Rs*X1thc + X1Rc*X1ths+ X1ths*dX1thc + X1thc*dX1ths)*dZ0 + X1thc*X1ths*ddZ0)/2/R0**2 - \
    (X1thc*X1ths*dR0*dZ0)/R0**3 

Z2c = X2zc - (X1thc*dX1zc-X1ths*dX1zs+2*X2thc*dZ0)/2/R0 + \
    (X1Rc*X1thc - X1Rs*X1ths + X1thc*dX1thc - X1ths*dX1ths)*dZ0/2/R0**2 - \
    + (X1thc**2 - X1ths**2)*ddZ0/4/R0**2 - (X1thc**2 - X1ths**2)*dR0*dZ0/2/R0**3)

dpanici avatar May 24 '24 18:05 dpanici