From 612d58baab0eae972398958a400cdd03533840a7 Mon Sep 17 00:00:00 2001 From: TQ Hirsch Date: Tue, 31 Jan 2023 16:16:04 +0100 Subject: [PATCH] Converted from roc to curvature --- app/Main.hs | 24 ++++++++++++++---------- lib/Petzval/Merit.hs | 31 ++++++++++++++++++++++++++++++- lib/Petzval/Optics.hs | 10 ++++++++-- lib/Petzval/Optics/RTM.hs | 4 ++-- lib/Petzval/System.hs | 8 ++++---- lib/Petzval/Trace.hs | 9 +++++++-- petzval-hs.cabal | 23 +++-------------------- 7 files changed, 68 insertions(+), 41 deletions(-) diff --git a/app/Main.hs b/app/Main.hs index 56c1f8c..30dacaf 100644 --- a/app/Main.hs +++ b/app/Main.hs @@ -15,7 +15,7 @@ import Petzval.Calculations import Petzval.Optimization import Petzval.Types import System.Environment (getArgs) -import System.IO (hPrint, stderr) +import System.IO (hPrint, hPutStrLn, stderr) import Linear import Numeric.AD.Mode (Scalar, Mode) import qualified Numeric.LinearAlgebra as L @@ -34,19 +34,19 @@ n_ssk8 = SellemeierMat [ (1.44857867, 1.17965926e-01) system1 = [ Stop{_thickness = 0, _outsideRadius=5} - , Surface{_material = bk7, _thickness = 10, _roc = 100, _outsideRadius=10} - , Surface{_material = air, _thickness = 95, _roc = -100, _outsideRadius=10} + , Surface{_material = bk7, _thickness = 10, _curvature = 1/100, _outsideRadius=10} + , Surface{_material = air, _thickness = 95, _curvature = -1/100, _outsideRadius=10} ] system2 = [ - Surface{_material = s_sk16, _outsideRadius=11.5, _roc=42.98790, _thickness = 4} - , Surface{_material = air, _outsideRadius = 11.5, _roc = -248.07740, _thickness = 10.51018} - , Surface{_material = s_f4, _outsideRadius = 9.852, _roc = -38.21035, _thickness = 2.5} - , Surface{_material = air, _outsideRadius = 8.885, _roc = 43.95894, _thickness = 0} + Surface{_material = s_sk16, _outsideRadius=11.5, _curvature=1/42.98790, _thickness = 4} + , Surface{_material = air, _outsideRadius = 11.5, _curvature = -1/248.07740, _thickness = 10.51018} + , Surface{_material = s_f4, _outsideRadius = 9.852, _curvature = -1/38.21035, _thickness = 2.5} + , Surface{_material = air, _outsideRadius = 8.885, _curvature = 1/43.95894, _thickness = 0} , Stop{_outsideRadius = 8.6762522794, _thickness = 9.86946} - , Surface{_material=s_sk16, _outsideRadius = 11, _roc=656.66349, _thickness = 4.5} - , Surface{_material = air, _outsideRadius = 11, _roc = -33.50754, _thickness = 86.48643} + , Surface{_material=s_sk16, _outsideRadius = 11, _curvature=1/656.66349, _thickness = 4.5} + , Surface{_material = air, _outsideRadius = 11, _curvature = -1/33.50754, _thickness = 86.48643} ] where s_sk16 = constMat 1.620408 :: SellemeierMat s_f4 = constMat 1.616589 @@ -75,9 +75,11 @@ doOptimize steps = runWriterT $ optimizeDLS cfg vars merit system1 , DynMerit $ TW 0 1 $ SpotSize Nothing , DynMerit $ TW 0 1 $ SpotSize (Just $ TraceConditions 7 0.587618) , DynMerit $ TW 0 1 $ SpotSize (Just $ TraceConditions 15 0.587618) + , DynMerit $ TW 1 100 $ Vign (TraceConditions 15 wl) ] + wl = 0.5875618 vars :: VariableSet - vars = (ix 1 . roc) `adjoin` (ix 2 . roc) `adjoin` (ix 2.thickness) + vars = (ix 1 . curvature) `adjoin` (ix 2 . curvature) `adjoin` (ix 2.thickness) cfg = def{maxSteps = steps,eps1=0,eps2=0} @@ -90,8 +92,10 @@ main = in do [arg] <- getArgs let steps = read arg + hPrint stderr system1 print $ L.norm_2 (L.fromList [3.0 , 4.0] :: L.Vector Double) (result, tr) <- doOptimize steps + hPutStrLn stderr $ "Result BFL: " ++ (show . rearFocalLength . lensRTM . bake 0.5875618) result hPrint stderr result putStrLn "---" mapM_ (putStrLn . DL.intercalate "," . map show) tr diff --git a/lib/Petzval/Merit.hs b/lib/Petzval/Merit.hs index 5669d9c..6c92d2c 100644 --- a/lib/Petzval/Merit.hs +++ b/lib/Petzval/Merit.hs @@ -1,5 +1,17 @@ {- LANGUAGE ImpredicativeTypes #-} -module Petzval.Merit (BFL(..), SpotSize(..), DynMerit(..), TW(..), Merit(..), MeritFunction, TraceConditions(..), evalMerit, mkMerit) where +module Petzval.Merit + ( + -- * Merit functions + BFL(..) + , SpotSize(..) + , DynMerit(..) + , Vign(..) + , TW(..) + , Merit(..) + , MeritFunction + , TraceConditions(..) + , evalMerit + , mkMerit) where import Petzval.Calculations import Petzval.Types @@ -9,6 +21,7 @@ import Petzval.System import Control.Lens import Data.Either import Data.Default +import Data.List (genericLength) import Petzval.Optics.RTM import Linear import Data.Maybe @@ -75,6 +88,22 @@ spotSize' conditions ts = rmsSize spotSize :: Calcuable a => MeritPart a spotSize = defTraceConditions >>= spotSize' +-- | Vignetting at a field angle. +-- Computes the proportion of rays that make it through the system +data Vign = Vign TraceConditions + deriving (Show) +instance Merit Vign where + calc (Vign cond) ts = (/ genericLength rays) + . genericLength + . rights + . map fst + . map (raytrace system) + $ rays + where rays = map (createRay Nothing ep fa) $ hexapolarPattern 6 + system = systemAtWavelength (wavelength cond) ts + fa = auto $ fieldAngle cond + ep = entrancePupil system + -- | Compute the back focal length bfl :: Calcuable a => MeritPart a -- | Compute the back focal length at a specific wavelength diff --git a/lib/Petzval/Optics.hs b/lib/Petzval/Optics.hs index 6b8681a..de1154a 100644 --- a/lib/Petzval/Optics.hs +++ b/lib/Petzval/Optics.hs @@ -16,6 +16,7 @@ module Petzval.Optics , outsideRadius , material , roc + , curvature , liftFp , specialize, bake ) where @@ -59,7 +60,7 @@ data Element mat fp = -- | Refractive surface Surface { _thickness :: fp , _outsideRadius :: fp - , _roc :: fp + , _curvature :: fp , _material :: mat } -- | Aperture stop @@ -67,9 +68,14 @@ data Element mat fp = _thickness :: fp , _outsideRadius :: fp } + -- | Imaging plane. This is a hack to make sure that rays go all the way to the target plane + | ImagingPlane { _thickness :: fp } deriving (Show) makeLenses ''Element +roc :: Fractional fp => Traversal' (Element mat fp) fp +roc = curvature . iso (1/) (1/) + -- | Determine if an element is a stop isStop :: Element mat a -> Bool isStop Stop{} = True @@ -82,7 +88,7 @@ isSurface _ = False -- | Translate a lens element from one FP type to another. E.g., can be used to convert from scalars to the types in an automatic differentiation tool. liftFp :: (Applicative f) => (fp -> f fp') -> Element m fp -> f (Element m fp') liftFp inj (s@Stop{_thickness, _outsideRadius}) = Stop <$> inj _thickness <*> inj _outsideRadius -liftFp inj (s@Surface{_thickness, _outsideRadius, _roc, _material}) = Surface <$> inj _thickness <*> inj _outsideRadius <*> inj _roc <*> pure _material +liftFp inj (s@Surface{_thickness, _outsideRadius, _curvature, _material}) = Surface <$> inj _thickness <*> inj _outsideRadius <*> inj _curvature <*> pure _material {-# INLINE liftFp #-} -- | Specialize a system for a specific wavelength diff --git a/lib/Petzval/Optics/RTM.hs b/lib/Petzval/Optics/RTM.hs index 391fdda..86b259a 100644 --- a/lib/Petzval/Optics/RTM.hs +++ b/lib/Petzval/Optics/RTM.hs @@ -40,10 +40,10 @@ systemRTM els = foldl1 (flip (!*!)) . map elementRTM $ els -- | Compute the refractive part of the RTM for an element refractionRTM :: (Fractional a, Mode a, Scalar a ~ Double) => Element BakedIOR a -> M22 a -refractionRTM (Surface{_thickness, _outsideRadius, _roc, _material=BakedIOR n1 n2}) = +refractionRTM (Surface{_thickness, _outsideRadius, _curvature, _material=BakedIOR n1 n2}) = let ior' = auto (n1 / n2) in V2 (V2 1 0) - (V2 ((ior' - fromIntegral 1) / _roc) ior') + (V2 ((ior' - fromIntegral 1) * _curvature) ior') refractionRTM Stop{} = V2 (V2 1 0) (V2 0 1) -- | Compute the part of the RTM that represents the thickness of the element diff --git a/lib/Petzval/System.hs b/lib/Petzval/System.hs index 504109e..6fa1488 100644 --- a/lib/Petzval/System.hs +++ b/lib/Petzval/System.hs @@ -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 = diff --git a/lib/Petzval/Trace.hs b/lib/Petzval/Trace.hs index 8ab6efa..5cc0c28 100644 --- a/lib/Petzval/Trace.hs +++ b/lib/Petzval/Trace.hs @@ -85,9 +85,10 @@ hitTest Stop{_outsideRadius} (Ray pos dir) = where dz = pos ^. _z / dir ^. _z npos = pos ^-^ (dir ^* dz) pass = quadrance (npos ^. _xy) <= _outsideRadius ^ 2 -hitTest Surface{_roc, _outsideRadius} ray@(Ray pos dir) = +hitTest Surface{_curvature, _outsideRadius} ray@(Ray pos dir) = toMaybe (hit1 && hit2) (Ray npos dir, Just normal) - where origin = pos & _z -~ _roc + where _roc = 1/_curvature + origin = pos & _z -~ _roc !a = dir `dot` dir !b = (dir `dot` origin) * 2 !c = (origin `dot` origin) - _roc ^ 2 @@ -103,6 +104,10 @@ hitTest Surface{_roc, _outsideRadius} ray@(Ray pos dir) = !normal = if (normal0 ^. _z < 0) then -normal0 else normal0 !npos = pos ^+^ dir ^* dist !hit2 = (quadrance $ npos ^. _xy) <= _outsideRadius^2 +hitTest ImagingPlane{_thickness} (Ray pos dir) = + Just (Ray npos dir, Nothing) + where dz = pos ^. _z / dir ^. _z - _thickness + npos = pos ^-^ (dir ^* dz) refract :: (Floating a, Ord a, Mode a, Scalar a ~ Double, Epsilon a) => BakedIOR -> V3 a -> Ray a -> Maybe (Ray a) refract (BakedIOR n1 n2) normal (Ray pos incident) = diff --git a/petzval-hs.cabal b/petzval-hs.cabal index f1012ef..bd5effe 100644 --- a/petzval-hs.cabal +++ b/petzval-hs.cabal @@ -44,25 +44,6 @@ executable petzval-hs data-default, petzval-hs hs-source-dirs: app - -executable petzval-prof - import: base - main-is: Main.hs - build-depends: base, - ad, - mtl, - linear, - lens, - -- criterion, - hmatrix, - data-default, - petzval-hs - hs-source-dirs: app - ghc-options: - -O2 - -threaded - -fprof-auto - -rtsopts library import: base @@ -78,6 +59,8 @@ library other-modules: Petzval.Internal.Vec hs-source-dirs: lib + ghc-options: + -O2 build-depends: base, lens, ad, @@ -88,4 +71,4 @@ library containers, monad-loops, hmatrix, - data-default + data-default \ No newline at end of file