Added README; added step to install BLAS in the build

.drone.yml

@@ -0,0 +1,11 @@
+kind: pipeline
+type: docker
+name: default
+
+steps:
+  - name: build
+    image: fpco/stack-build:lts-20
+    commands:
+      - apt-get update
+      - apt-get install -y libopenblas-dev
+      - stack build

README.md

@@ -0,0 +1,10 @@
+Petzval
+=======
+
+[![Build Status](https://ci.thequux.com/api/badges/thequux/petzval/status.svg)](https://ci.thequux.com/thequux/petzval)
+
+Petzval is a WIP lens optimization kernel. At this stage of development, it's
+probably only interesting to its author and people interested in contributing, 
+but we aim to provide an engine that can power something like Zemax or Code-V.
+
+

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

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

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

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

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 = 

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) =

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