commit 11cd3ee503da1b074618f39ab9940f290540a6f3 Author: TQ Hirsch Date: Fri Oct 7 18:59:44 2022 +0200 Initial commit diff --git a/.envrc b/.envrc new file mode 100644 index 0000000..3550a30 --- /dev/null +++ b/.envrc @@ -0,0 +1 @@ +use flake diff --git a/app/Main.hs b/app/Main.hs new file mode 100644 index 0000000..2385f93 --- /dev/null +++ b/app/Main.hs @@ -0,0 +1,45 @@ +module Main where + +import Petzval.Optics +import Petzval.Optics.RTM +import Petzval.System + +bk7 = SellemeierMat [ (1.03961212, 6.00069867e-3 ) + , (0.231792344, 2.00179144e-2 ) + , (1.01046945, 103.560653 )] +n_sk16 = SellemeierMat [ (1.343177740, 0.007046873) + , (0.241144399, 0.0229005000) + , (0.994317969, 92.75085260)] + +n_ssk8 = SellemeierMat [ (1.44857867, 1.17965926e-01) + , (1.06937528, 8.69310149E-03) + , (4.21566593E-02, 1.11300666E+02) ] + +system1 = + [ Stop{_thickness = 0, _outsideRadius=5} + , Surface{_material = bk7, _thickness = 10, _roc = 100, _outsideRadius=0} + , Surface{_material = air, _thickness = 95, _roc = -100, _outsideRadius=0} + ] + +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} + , 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} + ] + where s_sk16 = ConstMat 1.620408 + s_f4 = ConstMat 1.616589 + +main :: IO () +main = + let baked = bake 0.587618 system2 + ep = (entrancePupil baked){position=20.4747094540} + fa = 20 + in do + putStrLn $ show ep + putStrLn $ show $ rearFocalLength $ lensRTM baked + putStrLn $ show $ mconcat (seidel ep fa baked :: [Seidel Double]) diff --git a/flake.lock b/flake.lock new file mode 100644 index 0000000..1f429f4 --- /dev/null +++ b/flake.lock @@ -0,0 +1,41 @@ +{ + "nodes": { + "flake-utils": { + "locked": { + "lastModified": 1659877975, + "narHash": "sha256-zllb8aq3YO3h8B/U0/J1WBgAL8EX5yWf5pMj3G0NAmc=", + "owner": "numtide", + "repo": "flake-utils", + "rev": "c0e246b9b83f637f4681389ecabcb2681b4f3af0", + "type": "github" + }, + "original": { + "owner": "numtide", + "repo": "flake-utils", + "type": "github" + } + }, + "nixpkgs": { + "locked": { + "lastModified": 1661704917, + "narHash": "sha256-h1deRhxLw9iaYzIovHT9mLFuZq/9hoge+pXSbX98B78=", + "owner": "NixOS", + "repo": "nixpkgs", + "rev": "adf66cc31390f920854340679d7505ac577a5a8b", + "type": "github" + }, + "original": { + "id": "nixpkgs", + "type": "indirect" + } + }, + "root": { + "inputs": { + "flake-utils": "flake-utils", + "nixpkgs": "nixpkgs" + } + } + }, + "root": "root", + "version": 7 +} diff --git a/flake.nix b/flake.nix new file mode 100644 index 0000000..ae6b274 --- /dev/null +++ b/flake.nix @@ -0,0 +1,26 @@ +{ + description = "Flake utils demo"; + + inputs.flake-utils.url = "github:numtide/flake-utils"; + + outputs = { self, nixpkgs, flake-utils }: + flake-utils.lib.eachDefaultSystem (system: + let pkgs = nixpkgs.legacyPackages.${system}; in + { + packages = rec { + hello = pkgs.hello; + default = hello; + }; + apps = rec { + hello = flake-utils.lib.mkApp { drv = self.packages.${system}.hello; }; + default = hello; + }; + devShell = pkgs.mkShell { + buildInputs = with pkgs; [ + cabal-install + ghc + ]; + }; + } + ); +} diff --git a/lib/Petzval/Optics.hs b/lib/Petzval/Optics.hs new file mode 100644 index 0000000..a0e1b29 --- /dev/null +++ b/lib/Petzval/Optics.hs @@ -0,0 +1,105 @@ +{-# LANGUAGE TemplateHaskell #-} +{-# LANGUAGE NamedFieldPuns #-} +{-| +Basical optical types +-} +module Petzval.Optics + ( + -- * IOR calculations + Material(..), SellemeierMat(..) + , BakedIOR(..) + , ConstMat(..) + -- * Optical elements + , Element(..) + , thickness + , outsideRadius + , material + , roc + , liftFp + , specialize, bake + ) where + +import Control.Lens +import Data.Maybe + +-- | An optical material, suitable to calculate indices of refraction +class Material mat where + iorAtWavelength :: Double -> mat -> Double + -- | A material with a constant IOR of 1 + air :: mat + +-- | Calculates IOR based on the Sellemeier equation +data SellemeierMat = SellemeierMat [(Double, Double)] + deriving (Show) +instance Material SellemeierMat where + iorAtWavelength λ (SellemeierMat coeff) = sqrt $ 1 + (sum . map contrib $ coeff) + where contrib (b,c) = b * w2 / (w2 - c) + w2 = λ ^ 2 + air = SellemeierMat [] + +-- | A material with a constant IOR +newtype ConstMat = ConstMat Double + deriving (Eq, Ord, Show) + +instance Material ConstMat where + iorAtWavelength _ (ConstMat ior) = ior + air = ConstMat 1 + +-- | A representation of IOR as (previous,new) IOR at a surface +data BakedIOR = BakedIOR Double Double + deriving (Eq, Ord, Show) + +-- | An optical element +data Element mat fp = + -- | Refractive surface + Surface { _thickness :: fp + , _outsideRadius :: fp + , _roc :: fp + , _material :: mat + } + -- | Aperture stop + | Stop { + _thickness :: fp + , _outsideRadius :: fp + } + deriving (Show) +makeLenses ''Element + {- +-- | The space after the current element +thickness :: Lens (Element mat a) (Element mat a) a a +thickness inj (s@Surface{_thickness}) = (\nt -> s{_thickness=nt}) <$> inj _thickness +thickness inj (s@Stop{_thickness}) = (\nt -> s{_thickness=nt}) <$> inj _thickness +-- | The outside radius of the element. Rays that intersect the element beyond this radius are considered to have missed. +outsideRadius :: Lens (Element mat a) (Element mat a) a a +outsideRadius inj (s@Surface{_outsideRadius}) = (\nt -> s{_outsideRadius=nt}) <$> inj _outsideRadius +outsideRadius inj (s@Stop{_outsideRadius}) = (\nt -> s{_outsideRadius=nt}) <$> inj _outsideRadius + +infinity :: (RealFloat a) => a +infinity = encodeFloat (floatRadix 0 - 1) (snd $ floatRange 0) +-- | The radius of curvature of an element. This is 1 for stops +roc :: RealFloat a => Lens (Element mat a) (Element mat a) a a +roc inj (s@Surface{_roc}) = (\nt -> s{_roc=nt}) <$> inj _roc +roc inj (s@Stop{}) = (const s) <$> inj infinity + +-- | The material that a surface transitions into +-- +-- For stops, this returns None and ignores being set. +material :: (Material mat, Material mat') => Lens (Element mat a) (Element mat' a) (Maybe mat) (Maybe mat') +material inj (s@Surface{_thickness, _outsideRadius, _roc, _material}) = (\nt -> Surface{_material=fromMaybe air nt, _thickness, _outsideRadius, _roc}) <$> inj (Just _material) +material inj (s@Stop{_thickness, _outsideRadius}) = (const Stop{_thickness, _outsideRadius}) <$> inj Nothing +-} +-- | 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@Surface{_thickness, _outsideRadius, _roc, _material}) = (\t' or' roc' -> Surface{_thickness=t', _outsideRadius=or', _roc=roc', _material}) <$> inj _thickness <*> inj _outsideRadius <*> inj _roc +liftFp inj (s@Stop{_thickness, _outsideRadius}) = (\t' or' -> Stop{_thickness=t', _outsideRadius=or'}) <$> inj _thickness <*> inj _outsideRadius +{-# INLINE liftFp #-} + +specialize :: (Material mat) => Double -> [Element mat a] -> [Element Double a] +specialize wavelength = (each.material) %~ (iorAtWavelength wavelength) +bake :: (Material mat) => Double -> [Element mat a] -> [Element BakedIOR a] +bake wavelength mat = + snd $ + mapAccumLOf (each.material) (\n1 n2 -> (n2, BakedIOR n1 n2)) 1 $ + specialize wavelength mat + + diff --git a/lib/Petzval/Optics/RTM.hs b/lib/Petzval/Optics/RTM.hs new file mode 100644 index 0000000..391fdda --- /dev/null +++ b/lib/Petzval/Optics/RTM.hs @@ -0,0 +1,79 @@ +{-# LANGUAGE NamedFieldPuns #-} +{-| +Description: First-order optical computation tools based on ray transfer matrices + +Provides tools for computing various first-order properties of a lens +based on the paraxial approximation. +-} +module Petzval.Optics.RTM + ( + -- * Computation of ray transfer matrices from elements + elementRTM, refractionRTM, thicknessRTM + -- * Computation of ray transfer matrices for a system + , lensRTM, systemRTM + -- * First order properties of a lens system + , rearPrincipalPlane, rearFocalLength, rearWorkingDistance + -- * RTM manipulation + , reverseRTM + ) where + + +import Numeric.AD.Mode +import Linear.V2 +import Linear.Matrix +import Petzval.Optics +import Control.Lens + +-- | Compute an RTM for the optically active parts of a lens +-- system. This is the composition of the RTM for each element, but +-- ignores the thickness of the final element. +lensRTM :: (Fractional a, Mode a, Scalar a ~ Double) => [Element BakedIOR a] -> M22 a +-- | Compute the RTM of a complete optical system, including the +-- thickness of the final element. +systemRTM :: (Fractional a, Mode a, Scalar a ~ Double) => [Element BakedIOR a] -> M22 a +lensRTM [] = V2 (V2 1 0) (V2 0 1) +lensRTM [el] = refractionRTM el +lensRTM (el:rest) = lensRTM rest !*! elementRTM el +systemRTM [] = V2 (V2 1 0) (V2 0 1) +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}) = + let ior' = auto (n1 / n2) + in V2 (V2 1 0) + (V2 ((ior' - fromIntegral 1) / _roc) ior') +refractionRTM Stop{} = V2 (V2 1 0) (V2 0 1) + +-- | Compute the part of the RTM that represents the thickness of the element +thicknessRTM :: (Fractional a) => Element mat a -> M22 a +thicknessRTM element = V2 (V2 1 (element ^.thickness)) (V2 0 1) + +-- | Compute the complete RTM of an element, i.e., \(thickness * refractive\) +elementRTM :: (Fractional a, Mode a, Scalar a ~ Double) => Element BakedIOR a -> M22 a +elementRTM el@Surface{} = thicknessRTM el !*! refractionRTM el +elementRTM el@Stop{} = thicknessRTM el + +-- | Compute the position of the rear principal plane for a lens system relative to the final element. +rearPrincipalPlane :: Fractional a => M22 a -> a +rearPrincipalPlane rtm = let V2 y m = rtm !* V2 1 0 + in (1 - y) / m +-- | Compute the rear focal length (distance from the rear principal plane to the focal point) +rearFocalLength :: Fractional a => M22 a -> a +rearFocalLength rtm = let V2 y m = rtm !* V2 1 0 + in -1 / m +-- | Compute the rear working distance of a system (distance from the final element to the focal point) +rearWorkingDistance :: Fractional a => M22 a -> a +rearWorkingDistance rtm = let V2 y m = rtm !* V2 1 0 + in -y / m + +-- | Compute the RTM for the reverse of the system represented by @rtm@. This is equivalent to +-- \[ +-- \begin{bmatrix} 1 & 0 \\ 0 & -1 \end{bmatrix} +-- rtm^{-1} +-- \begin{bmatrix} 1 & 0 \\ 0 & -1 \end{bmatrix} +-- \] +reverseRTM :: (Fractional a) => M22 a -> M22 a +reverseRTM rtm = let V2 (V2 a b) (V2 c d) = inv22 rtm + in V2 (V2 a (-b)) (V2 (-c) d) diff --git a/lib/Petzval/System.hs b/lib/Petzval/System.hs new file mode 100644 index 0000000..99ac9f9 --- /dev/null +++ b/lib/Petzval/System.hs @@ -0,0 +1,126 @@ +module Petzval.System + ( Pupil(..) + , entrancePupil + , createRay + , seidel, Seidel(..) + , Ray + ) where + +import Petzval.Optics +import Petzval.Optics.RTM +import Linear +import Numeric.AD.Mode +import Control.Lens +import Data.List + +import qualified Debug.Trace + +data Ray a = Ray (V3 a) (V3 a) + +splitSystem :: [ Element mat a ] -> ([Element mat a], Maybe (Element mat a), [Element mat a]) +splitSystem (s@Stop{}:rest) = ([], Just s, rest) +splitSystem [] = ([], Nothing, []) +splitSystem (s:rest) = (s:pfx, stop, sfx) + where (pfx, stop, sfx) = splitSystem rest + +data Pupil a = Pupil { radius :: a, position :: a } + deriving (Show) + +pupil (stop@Stop{_outsideRadius}) subsystem = + let rtm = systemRTM subsystem + V2 my mu = rtm !* V2 1 0 + V2 cy cu = (inv22 rtm) !* V2 0 (-1) + --msg = intercalate " " ["ptrace: ", show (rtm !* V2 (_outsideRadius / my ) 0)] + in Pupil { radius = _outsideRadius / my + , position = -cy / cu } +entrancePupil :: (RealFloat a, Mode a, Scalar a ~ Double, Epsilon a, Show a) => [Element BakedIOR a] -> Pupil a +entrancePupil system = pupil stop prefix + where (prefix, (Just stop), _) = splitSystem system + +createRay :: (RealFloat a, Mode a, Scalar a ~ Double, Epsilon a) => Maybe a -> Pupil a -> a -> (a,a) -> Ray a +createRay (Just objectPlane) Pupil{position=pz,radius=pr} h (px, py) = + Ray source (normalize $ target ^-^ source) + where dz = pz - objectPlane + source = V3 0 (dz * tan h) objectPlane + target = V3 (px * pr) (py * pr) pz +createRay Nothing Pupil{position=pz,radius=pr} h (px, py) = Ray source (normalize $ target ^-^ source) + where h' = (pi * (-abs h) / 180) -- field angle in rad + dy = (V3 0 (cos h') (-sin h')) `project` (V3 0 (py * pr) 0) + dz = V3 0 (pz * tan h') (pz * cos h') + + source = dy ^+^ dz + target = V3 (px * pr) (py * pr) pz + +data Seidel a = Seidel + { sphr, coma, asti, fcur, dist :: a + -- , c1, c2 :: a + } + deriving (Show) + +instance Num a => Semigroup (Seidel a) where + (<>) Seidel{sphr, coma, asti, fcur, dist} Seidel{sphr=sphr', coma=coma', asti=asti', fcur=fcur', dist=dist'} = + Seidel { sphr = sphr + sphr' + , coma = coma + coma' + , asti = asti + asti' + , fcur = fcur + fcur' + , dist = dist + dist' + --, c1 = c1 + c1' + --, c2 = c2 + c2' + } + +instance Num a => Monoid (Seidel a) where + mempty = Seidel { sphr=0, coma=0, asti=0, fcur=0, dist=0 + -- , c1 = mempty, c2 = mempty + } + mappend Seidel{sphr, coma, asti, fcur, dist} Seidel{sphr=sphr', coma=coma', asti=asti', fcur=fcur', dist=dist'} = + Seidel { sphr = sphr + sphr' + , coma = coma + coma' + , asti = asti + asti' + , fcur = fcur + fcur' + , dist = dist + dist' + --, c1 = c1 + c1' + --, c2 = c2 + c2' + } + + + +-- | 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' 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}) + where rays' = refractionRTM s !*! rays + rays'' = thicknessRTM s !*! rays' + marg = column _y + chief = column _x + _i = to (\ray -> (ray ^. _x) / roc + (ray ^. _y)) + a = auto n1 * (rays ^. marg._i) + abar = auto n1 * (rays ^. chief._i) + h = rays ^. marg._x + δun = (rays' ^. marg._y / auto n2) - (rays ^. marg._y / auto n1) + δ1n = auto (1/n2-1/n1) + δ1n2 = auto (1/n2^2-1/n1^2) + lagrange = auto n1 * det22 rays + + sphr = a ^ 2 * h * δun + coma = -a * abar * h * δun + asti = -abar^2 * h * δun + fcur = lagrange^2 / roc * δ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 + msg = intercalate "\t" [show δun, show a, show abar] + -- TODO: evaluate at other IORs + -- cbas = h (n2 - n1) / n1 + -- c1 = + + + + +seidel :: (RealFloat a, Mode a, Scalar a ~ Double, Epsilon a, Show a) => Pupil a -> Double -> [Element BakedIOR a] -> [Seidel a] +seidel pupil fieldAngle = + snd . mapAccumL seidel' rays + where ubar = auto (tan (pi / 180 * fieldAngle)) + rays = V2 (V2 (-position pupil * ubar) (radius pupil)) + (V2 ubar 0) + diff --git a/petzval-hs.cabal b/petzval-hs.cabal new file mode 100644 index 0000000..9ed8736 --- /dev/null +++ b/petzval-hs.cabal @@ -0,0 +1,53 @@ +cabal-version: 2.4 +name: petzval-hs +version: 0.1.0.0 + +-- A short (one-line) description of the package. +-- synopsis: + +-- A longer description of the package. +-- description: + +-- A URL where users can report bugs. +-- bug-reports: + +-- The license under which the package is released. +-- license: +author: TQ Hirsch +maintainer: thequux@thequux.com + +-- A copyright notice. +-- copyright: +-- category: +extra-source-files: CHANGELOG.md + +common base + default-extensions: GADTs, NamedFieldPuns + default-language: Haskell2010 + +executable petzval-hs + import: base + main-is: Main.hs + + -- Modules included in this executable, other than Main. + -- other-modules: + + -- LANGUAGE extensions used by modules in this package. + -- other-extensions: + build-depends: base ^>=4.15.1.0, + petzval-hs + hs-source-dirs: app + + +library + import: base + exposed-modules: + Petzval.Optics + Petzval.Optics.RTM + Petzval.System + hs-source-dirs: lib + build-depends: base ^>=4.15.1.0, + lens ^>=5.2, + ad ^>=4.5.2, + linear ^>=1.21, + \ No newline at end of file