--
-- LUA script to compute multipoles in FEMM
--
-- Few standard cases are considered:
--  * dipole 180 deg (ex. C shape)
--  * dipole 90 deg (ex. H shape)
--  * quadrupole 45 deg (ex. standard symmetric quadrupole)
--
-- In all cases, the center is 0, 0 and the skew coefficients are 0 
-- 
-- The script computes two sets of multipoles:
--  * one from A (the vector potential)
--  * another one from a radial projection of B
-- They should be the same, so the difference in a way shows 
-- how much to trust these numbers; the ones from A should be better, 
-- as this is the finite element solution without further manipulations
-- (derivation, radial projection) while B is rougher (linear elements, 
-- so B is constant over each triangle), but then it's smoothed out in the postprocessor
--

case_index = 1
-- 1 ===> dipole 180 deg (ex. C shape)
-- 2 ===> dipole 90 deg (ex. H shape)
-- 3 ===> quadrupole 45 deg (ex. standard symmetric quadrupole)

nh = 15 -- number of harmonics
np = 256 -- number of samples points (remember Nyquist...)
R = 2*25/3 -- reference radius
Rs = 0.95*25 -- sampling radius, can be the same as R or the largest still in the air

if case_index == 1 then
   thmax = pi
   ihmin, ihstep, ihfund = 1, 1, 1
elseif case_index == 2 then
   thmax = pi/2
   ihmin, ihstep, ihfund = 1, 2, 1
elseif case_index == 3 then
   thmax = pi/4
   ihmin, ihstep, ihfund = 2, 4, 2
end

-- multiplication factor, to include mirror copies
fact = 2*pi/thmax 

-- angular increment
dth = thmax/(np-1)

-- sampling Br and A
Br, Ath = {}, {}
A, B1, B2 = mo_getpointvalues(0, 0)
Actr = A
for ip = 1, np do
   th = (ip-1)*dth
   xs, ys = Rs*cos(th), Rs*sin(th)
   A, B1, B2 = mo_getpointvalues(xs, ys)
   Br[ip] = B1*cos(th) + B2*sin(th)
   Ath[ip] = A - Actr
end

-- print header, with reference and sampling radii
print('Reference radius ', R)
print('Sampling radius ', Rs)

-- harmonics from Br
B, b = {}, {}
for ih = 1, nh do
   B[ih] = 0
end
for ih = ihmin, nh, ihstep do
   B[ih] = Br[1]*sin(0)/2
   for ip = 2, np-1 do
      B[ih] = B[ih] + Br[ip]*sin(ih*((ip-1)*dth))
   end
   B[ih] = B[ih] + Br[np]*sin(ih*thmax)/2
   B[ih] = fact*B[ih]*dth/pi
   B[ih] = (R/Rs)^(ih-1)*B[ih]
end
print('from Br')
print(B[ihfund])
for ih = ihmin, nh, ihstep do
   b[ih] = B[ih]/B[ihfund]*10000
   print(ih, b[ih])
end

-- harmonics from A
B, b = {}, {}
for ih = 1, nh do
   B[ih] = 0
end
for ih = ihmin, nh, ihstep do
   B[ih] = Ath[1]*cos(0)/2
   for ip = 2, np-1 do
      B[ih] = B[ih] + Ath[ip]*cos(ih*((ip-1)*dth))
   end
   B[ih] = B[ih] + Ath[np]*cos(ih*thmax)/2
   B[ih] = fact*B[ih]*dth/pi*(-ih/Rs)*1000
   B[ih] = (R/Rs)^(ih-1)*B[ih]
end
print('from A')
print(B[ihfund])
for ih = ihmin, nh, ihstep do
   b[ih] = B[ih]/B[ihfund]*10000
   print(ih, b[ih])
end