Converted from roc to curvature

This commit is contained in:
2023-01-31 16:16:04 +01:00
parent 0bd2ef581c
commit 612d58baab
7 changed files with 68 additions and 41 deletions

View File

@@ -65,12 +65,12 @@ instance Num a => Monoid (Seidel a) where
-- | Initial matrix is [ h_ h; u_ u ]
seidel' :: (RealFloat a, Mode a, Scalar a ~ Double, Epsilon a) => M22 a -> Element BakedIOR a -> (M22 a, Seidel a)
seidel' rays (s@Stop{}) = (thicknessRTM s !*! rays, mempty)
seidel' rays (s@Surface{_material=BakedIOR n1 n2,_roc=roc}) = (rays'', Seidel {sphr,coma,asti,fcur,dist})
seidel' rays (s@Surface{_material=BakedIOR n1 n2,_curvature=curvature}) = (rays'', Seidel {sphr,coma,asti,fcur,dist})
where rays' = refractionRTM s !*! rays
rays'' = thicknessRTM s !*! rays'
marg = column _y
chief = column _x
_i = to (\ray -> (ray ^. _x) / roc + (ray ^. _y))
_i = to (\ray -> (ray ^. _x) / curvature + (ray ^. _y))
a = auto n1 * (rays ^. marg._i)
abar = auto n1 * (rays ^. chief._i)
h = rays ^. marg._x
@@ -82,11 +82,11 @@ seidel' rays (s@Surface{_material=BakedIOR n1 n2,_roc=roc}) = (rays'', Seidel {s
sphr = a ^ 2 * h * δun
coma = -a * abar * h * δun
asti = -abar^2 * h * δun
fcur = lagrange^2 / roc * δ1n
fcur = lagrange^2 * curvature * δ1n
-- Normally this is defined as abar/a * (asti*fcur)
-- However, a is 0 whenever the marginal ray is normal to the surface
-- Thus, we use this formula from Kidger
dist = -abar^3 * h * δ1n2 + (rays^. chief._x) * abar * (2 * h * abar - rays ^. chief._x * a) * δ1n / roc
dist = -abar^3 * h * δ1n2 + (rays^. chief._x) * abar * (2 * h * abar - rays ^. chief._x * a) * δ1n * curvature
-- TODO: evaluate at other IORs
-- cbas = h (n2 - n1) / n1
-- c1 =