First pass at merit function

This commit is contained in:
2022-10-10 11:46:27 +02:00
parent a80a7f749e
commit a61f15df44
5 changed files with 32 additions and 6 deletions

View File

@@ -1,8 +1,14 @@
module Main where module Main where
import Control.Lens
import Data.Either
import Petzval.Optics import Petzval.Optics
import Petzval.Optics.RTM import Petzval.Optics.RTM
import Petzval.System import Petzval.System
import Petzval.Trace
import Petzval.Calculations
import Linear
import Numeric.AD.Mode (Scalar, Mode)
bk7 = SellemeierMat [ (1.03961212, 6.00069867e-3 ) bk7 = SellemeierMat [ (1.03961212, 6.00069867e-3 )
, (0.231792344, 2.00179144e-2 ) , (0.231792344, 2.00179144e-2 )
@@ -34,6 +40,22 @@ system2 =
where s_sk16 = ConstMat 1.620408 where s_sk16 = ConstMat 1.620408
s_f4 = ConstMat 1.616589 s_f4 = ConstMat 1.616589
merit :: (Material mat, RealFloat a, Scalar a ~ Double, Epsilon a, Mode a) => [Element mat a] -> a
-- merit :: (Material mat) => [Element mat Double] -> Double
merit system = result
where baked = bake 0.587618 system
ep = entrancePupil baked
result = sum
. map (^2)
. map (\ angle -> rmsSize
. toListOf (each._pos._xy)
-- . (\x -> x :: [Ray a])
. rights
. map fst
. map (raytrace baked . createRay Nothing ep angle)
$ spiralPattern 100)
$ [0, 7.07, 15]
main :: IO () main :: IO ()
main = main =
let baked = bake 0.587618 system2 let baked = bake 0.587618 system2

View File

@@ -18,6 +18,7 @@
devShell = pkgs.mkShell { devShell = pkgs.mkShell {
buildInputs = with pkgs; [ buildInputs = with pkgs; [
cabal-install cabal-install
haskell-language-server
ghc ghc
]; ];
}; };

View File

@@ -25,7 +25,7 @@ class VConcat (v1 :: * -> *) (v2 :: * -> *) where
type VConcated v1 v2 :: * -> * type VConcated v1 v2 :: * -> *
vconcat :: v1 a -> v2 a -> (VConcated v1 v2) a vconcat :: v1 a -> v2 a -> (VConcated v1 v2) a
instance {-# OVERLAPPED #-} VConcat V0 v where instance VConcat V0 v where
type VConcated V0 v = v type VConcated V0 v = v
vconcat V0 v = v vconcat V0 v = v

View File

@@ -34,7 +34,7 @@ pupil (stop@Stop{_outsideRadius}) subsystem =
, position = -cy / cu } , position = -cy / cu }
-- | Compute the entrance pupil of a system viewed from infinity -- | Compute the entrance pupil of a system viewed from infinity
entrancePupil :: (RealFloat a, Mode a, Scalar a ~ Double, Epsilon a, Show a) => [Element BakedIOR a] -> Pupil a entrancePupil :: (RealFloat a, Mode a, Scalar a ~ Double, Epsilon a) => [Element BakedIOR a] -> Pupil a
entrancePupil system = pupil stop prefix entrancePupil system = pupil stop prefix
where (prefix, (Just stop), _) = splitSystem system where (prefix, (Just stop), _) = splitSystem system
@@ -73,9 +73,9 @@ instance Num a => Monoid (Seidel a) where
-- | Initial matrix is [ h_ h; u_ u ] -- | Initial matrix is [ h_ h; u_ u ]
seidel' :: (RealFloat a, Mode a, Scalar a ~ Double, Epsilon a, Show a) => M22 a -> Element BakedIOR a -> (M22 a, Seidel a) 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@Stop{}) = (thicknessRTM s !*! rays, mempty)
seidel' rays (s@Surface{_material=BakedIOR n1 n2,_roc=roc}) = Debug.Trace.trace msg (rays'', Seidel {sphr,coma,asti,fcur,dist}) seidel' rays (s@Surface{_material=BakedIOR n1 n2,_roc=roc}) = (rays'', Seidel {sphr,coma,asti,fcur,dist})
where rays' = refractionRTM s !*! rays where rays' = refractionRTM s !*! rays
rays'' = thicknessRTM s !*! rays' rays'' = thicknessRTM s !*! rays'
marg = column _y marg = column _y
@@ -97,13 +97,12 @@ seidel' rays (s@Surface{_material=BakedIOR n1 n2,_roc=roc}) = Debug.Trace.trace
-- However, a is 0 whenever the marginal ray is normal to the surface -- However, a is 0 whenever the marginal ray is normal to the surface
-- Thus, we use this formula from Kidger -- 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 / roc
msg = intercalate "\t" [show δun, show a, show abar]
-- TODO: evaluate at other IORs -- TODO: evaluate at other IORs
-- cbas = h (n2 - n1) / n1 -- cbas = h (n2 - n1) / n1
-- c1 = -- c1 =
-- | Compute the seidel coefficients of a system given the pupil and field angle -- | Compute the seidel coefficients of a system given the pupil and field angle
seidel :: (RealFloat a, Mode a, Scalar a ~ Double, Epsilon a, Show a) => Pupil a -> Double -> [Element BakedIOR a] -> [Seidel a] seidel :: (RealFloat a, Mode a, Scalar a ~ Double, Epsilon a) => Pupil a -> Double -> [Element BakedIOR a] -> [Seidel a]
seidel pupil fieldAngle = seidel pupil fieldAngle =
snd . mapAccumL seidel' rays snd . mapAccumL seidel' rays
where ubar = auto (tan (pi / 180 * fieldAngle)) where ubar = auto (tan (pi / 180 * fieldAngle))

View File

@@ -35,6 +35,9 @@ executable petzval-hs
-- LANGUAGE extensions used by modules in this package. -- LANGUAGE extensions used by modules in this package.
-- other-extensions: -- other-extensions:
build-depends: base ^>=4.15.1.0, build-depends: base ^>=4.15.1.0,
ad ^>=4.5.2,
linear ^>=1.21,
lens ^>=5.2,
petzval-hs petzval-hs
hs-source-dirs: app hs-source-dirs: app
@@ -46,6 +49,7 @@ library
Petzval.Optics.RTM Petzval.Optics.RTM
Petzval.System Petzval.System
Petzval.Trace Petzval.Trace
Petzval.Calculations
other-modules: other-modules:
Petzval.Internal.Vec Petzval.Internal.Vec
hs-source-dirs: lib hs-source-dirs: lib