diff --git a/0_hello_deeplens.py b/0_hello_deeplens.py
index eb981c2..2f3825e 100644
--- a/0_hello_deeplens.py
+++ b/0_hello_deeplens.py
@@ -17,10 +17,40 @@
 
 
 def main():
-    lens = GeoLens(filename="./lenses/camera/ef35mm_f2.0.json")
+    lens = GeoLens(filename="./lenses/camera/anamorphic_50mm.json")
     # lens = GeoLens(filename='./lenses/cellphone/cellphone80deg.json')
     # lens = GeoLens(filename='./lenses/zemax_double_gaussian.zmx')
-    lens.analysis(render=True)
+    lens.analysis(
+        f"./initial_lens",
+        zmx_format=True,
+        multi_plot=False,
+    )
+    lens.optimize(
+         lrs=[6e-4, 1e-4, 0.1, 1e-4],
+         decay=0.02,
+         iterations=5000,
+         centroid=True,
+         importance_sampling=True,
+         optim_mat=True,
+         shape_control=True,
+         anamorphic=True,
+         test_per_iter=20,
+         result_dir="./result",
+    )
+
+    # =====> 3. Analyze final result
+    lens.prune_surf(expand_surf=0.02)
+    lens.post_computation()
+
+    logging.info(
+        f"Actual: diagonal FOV {lens.hfov}, r sensor {lens.r_sensor}, F/{lens.fnum}."
+    )
+    lens.write_lens_json(f"{result_dir}/final_lens.json")
+    lens.analysis(save_name=f"{result_dir}/final_lens", zmx_format=True)
+
+    # =====> 4. Create video
+    create_video_from_images(f"{result_dir}", f"{result_dir}/autolens.mp4", fps=10)
+
 
 
 if __name__ == "__main__":
diff --git a/1_end2end_5lines.py b/1_end2end_5lines.py
index 367c000..127fe5e 100644
--- a/1_end2end_5lines.py
+++ b/1_end2end_5lines.py
@@ -66,11 +66,11 @@ def config():
         raise Exception("Add your wandb logging config here.")
 
     # ==> Device
-    num_gpus = torch.cuda.device_count()
+    num_gpus = 1
     args["num_gpus"] = num_gpus
-    device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
+    device = torch.device("cuda" if torch.cuda.is_available() else "mps")
     args["device"] = device
-    logging.info(f"Using {num_gpus} {torch.cuda.get_device_name(0)} GPU(s)")
+    #logging.info(f"Using {num_gpus} {torch.cuda.get_device_name(0)} GPU(s)")
 
     # ==> Save config and original code
     with open(f"{result_dir}/config.yml", "w") as f:
@@ -100,6 +100,7 @@ def end2end_train(lens, net, args):
 
     train_set = ImageDataset(args["train"]["train_dir"], lens.sensor_res)
     train_loader = DataLoader(train_set, batch_size=args["train"]["bs"])
+    logging.info(f'train_loader: {train_loader}')
 
     # ==> Network optimizer
     batchs = len(train_loader)
@@ -134,18 +135,22 @@ def end2end_train(lens, net, args):
         # ==> Train 1 epoch
         for img_org in tqdm(train_loader):
             img_org = img_org.to(device)
+            logging.info(f'img_org: {img_org}')
 
             # => Render image
             # ========================================
             # Line 3: plug-and-play diff-rendering
             # ========================================
             img_render = lens.render(img_org)
+            logging.info(f'img_render: {img_render}')
 
             # => Image restoration
             img_rec = net(img_render)
+            logging.info(f'img_rec: {img_rec}')
 
             # => Loss
             L_rec = cri_l1(img_rec, img_org)
+            logging.info(f'L_rec: {L_rec}')
 
             # => Back-propagation
             net_optim.zero_grad()
@@ -225,9 +230,11 @@ def end2end_train(lens, net, args):
     # Line 1: load a lens
     # ========================================
     lens = GeoLens(filename=args["lens"]["path"])
+    logging.info(f'lens: {lens}')
     lens.change_sensor_res(args["train"]["img_res"])
     net = UNet()
     net = net.to(lens.device)
+    logging.info(f'net: {net}')
     if args["network"]["pretrained"]:
         net.load_state_dict(torch.load(args["network"]["pretrained"]))
 
diff --git a/2_autolens_rms.py b/2_autolens_rms.py
index 850531a..b69181e 100644
--- a/2_autolens_rms.py
+++ b/2_autolens_rms.py
@@ -58,11 +58,11 @@ def config():
     logging.info(f"EXP: {args['EXP_NAME']}")
 
     # Device
-    num_gpus = torch.cuda.device_count()
+    num_gpus = 1 #torch.cuda.device_count()
     args["num_gpus"] = num_gpus
-    device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
+    device = torch.device("mps" if torch.backends.mps.is_available() else "cpu")
     args["device"] = device
-    logging.info(f"Using {num_gpus} {torch.cuda.get_device_name(0)} GPU(s)")
+    #logging.info(f"Using {num_gpus} {torch.cuda.get_device_name(0)} GPU(s)")
 
     # ==> Save config and original code
     with open(f"{result_dir}/config.yml", "w") as f:
@@ -153,9 +153,13 @@ def curriculum_design(
             center_p = center_p.unsqueeze(-2).repeat(1, 1, spp, 1)
 
         # =======================================
-        # Optimize lens by minimizing rms
+        # Optimize lens by minimizing rms and anamorphic squeeze
         # =======================================
         loss_rms = []
+        loss_anamorphic_list = []
+        delta = 1  # object-space offset (ensure consistent units)
+        target_squeeze = 1.5  # desired horizontal/vertical magnification ratio
+
         for j, wv in enumerate(WAVE_RGB):
             # Ray tracing to sensor
             ray = rays_backup[j].clone()
@@ -185,13 +189,40 @@ def curriculum_design(
             )
             loss_rms.append(l_rms)
 
+            # ---- Anamorphic (Squeeze) Loss Computation ----
+            # Horizontal offset ray.
+            ray_offset_h = self.get_cached_rays(depth=depth, wvln=wv, offset='h')
+            ray_offset_h = self.trace2sensor(ray_offset_h)
+            xy_offset_h = ray_offset_h.o[..., :2]
+            ra_xy_offset_h = ray_offset_h.ra.clone().detach()
+            xy_offset_h_norm = (xy_offset_h - center_p) * ra_xy_offset_h.unsqueeze(-1)
+
+            # Vertical offset ray.
+            ray_offset_v = self.get_cached_rays(depth=depth, wvln=wv, offset='v')
+            ray_offset_v = self.trace2sensor(ray_offset_v)
+            xy_offset_v = ray_offset_v.o[..., :2]
+            ra_xy_offset_v = ray_offset_v.ra.clone().detach()
+            xy_offset_v_norm = (xy_offset_v - center_p) * ra_xy_offset_v.unsqueeze(-1)
+
+            # Compute effective magnifications.
+            M_x = torch.mean(torch.abs(xy_offset_h_norm[..., 0])) / delta
+            M_y = torch.mean(torch.abs(xy_offset_v_norm[..., 1])) / delta
+            print(f'Magnifications: {M_x}, {M_y}')
+            ratio = M_x / (M_y + EPSILON)
+            loss_anamorphic_list.append((ratio - target_squeeze) ** 2)
+
         # RMS loss for all wavelengths
         loss_rms = sum(loss_rms) / len(loss_rms)
+        #loss_anamorphic = torch.mean(torch.stack(loss_anamorphic_list))
+        loss_anamorphic = sum(loss_anamorphic_list) / len(loss_anamorphic_list)
 
         # Lens design constraint
         loss_reg = self.loss_reg()
         w_reg = 0.1
-        L_total = loss_rms + w_reg * loss_reg
+        # Adding for anamorphics. 1.5x squeeze hardcoded for now
+        w_anamorphic = 1.0
+        print(f"Losses: {loss_rms}, {w_reg * loss_reg}, {w_anamorphic * loss_anamorphic}")
+        L_total = loss_rms + w_reg * loss_reg + w_anamorphic * loss_anamorphic
 
         # Gradient-based optimization
         optimizer.zero_grad()
@@ -236,7 +267,7 @@ def curriculum_design(
         lrs=[float(lr) for lr in args["lrs"]],
         decay=float(args["decay"]),
         iterations=5000,
-        test_per_iter=50,
+        test_per_iter=3,
         optim_mat=True,
         match_mat=False,
         shape_control=True,
diff --git a/6_hybridlens_design.py b/6_hybridlens_design.py
index 080bd5b..47f04b7 100644
--- a/6_hybridlens_design.py
+++ b/6_hybridlens_design.py
@@ -55,11 +55,11 @@ def config():
         raise Exception("Add your wandb logging config here.")
 
     # ==> Device
-    num_gpus = torch.cuda.device_count()
+    num_gpus = torch.cuda.device_count() if torch.cuda.is_available() else 0
     args["num_gpus"] = num_gpus
     device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
     args["device"] = device
-    logging.info(f"Using {num_gpus} {torch.cuda.get_device_name(0)} GPU(s)")
+    # logging.info(f"Using {num_gpus} {torch.cuda.get_device_name(0)} GPU(s)")
 
     # ==> Save config
     with open(f"{result_dir}/config.yml", "w") as f:
@@ -76,7 +76,7 @@ def main(args):
     # Create a hybrid refractive-diffractive lens
     lens = HybridLens(lens_path="./lenses/hybridlens/a489_doe.json")
     lens.refocus(foc_dist=-1000.0)
-    lens.double()
+    #lens.double()
 
     # PSF optimization loop to focus blue light
     optimizer = lens.get_optimizer(doe_lr=0.1, lens_lr=[1e-4, 1e-4, 1e-1, 1e-5])
diff --git a/configs/2_auto_lens_design.yml b/configs/2_auto_lens_design.yml
index 5b875aa..b50d0ae 100644
--- a/configs/2_auto_lens_design.yml
+++ b/configs/2_auto_lens_design.yml
@@ -2,16 +2,16 @@ DEBUG: True
 seed: ~
 
 # experiment settings
-EXP_NAME: 'Auto lens design'
+EXP_NAME: 'Auto lens design - anamorphic'
 
 # lens target example 1 (camera lens)
-foclen: 85.0
-fov: 40.0
+foclen: 50.0
+fov: 48.0
 fnum: 4.0
-flange: 18.0
-thickness: 120.0
-lens_type: [["Spheric", "Spheric"], ["Spheric", "Spheric"], ["Spheric", "Spheric", "Spheric"], ["Aperture"], ["Spheric", "Spheric"], ["Spheric", "Aspheric"], ["Spheric", "Aspheric"]]
-lrs: [5e-4, 1e-3, 1e-1, 1e-3]
+flange: 27.75
+thickness: 111.906
+lens_type: [["Anamorphic", "Anamorphic", "Anamorphic"], ["Spheric", "Spheric", "Spheric"], ["Spheric", "Spheric"], ["Aperture"], ["Spheric", "Spheric"], ["Spheric", "Spheric", "Spheric"], ["Anamorphic", "Anamorphic"], ["Anamorphic", "Anamorphic"]]
+lrs: [2e-3, 2e-3, 1e-1, 1e-3]
 decay: 0.001
 
 # # lens target example 2 (mobile lens)
@@ -22,4 +22,4 @@ decay: 0.001
 # thickness: 9.0
 # lens_type: [["Aperture"], ["Aspheric", "Aspheric"], ["Aspheric", "Aspheric"], ["Aspheric", "Aspheric"], ["Aspheric", "Aspheric"], ["Aspheric", "Aspheric"]]
 # lrs: [3e-4, 1e-4, 1e-1, 1e-2]
-# decay: 0.01
\ No newline at end of file
+# decay: 0.01
diff --git a/configs/2_auto_lens_design_original.yml b/configs/2_auto_lens_design_original.yml
new file mode 100644
index 0000000..5b875aa
--- /dev/null
+++ b/configs/2_auto_lens_design_original.yml
@@ -0,0 +1,25 @@
+DEBUG: True
+seed: ~
+
+# experiment settings
+EXP_NAME: 'Auto lens design'
+
+# lens target example 1 (camera lens)
+foclen: 85.0
+fov: 40.0
+fnum: 4.0
+flange: 18.0
+thickness: 120.0
+lens_type: [["Spheric", "Spheric"], ["Spheric", "Spheric"], ["Spheric", "Spheric", "Spheric"], ["Aperture"], ["Spheric", "Spheric"], ["Spheric", "Aspheric"], ["Spheric", "Aspheric"]]
+lrs: [5e-4, 1e-3, 1e-1, 1e-3]
+decay: 0.001
+
+# # lens target example 2 (mobile lens)
+# foclen: 6.0
+# fov: 70.0
+# fnum: 2.0
+# flange: 1.0
+# thickness: 9.0
+# lens_type: [["Aperture"], ["Aspheric", "Aspheric"], ["Aspheric", "Aspheric"], ["Aspheric", "Aspheric"], ["Aspheric", "Aspheric"], ["Aspheric", "Aspheric"]]
+# lrs: [3e-4, 1e-4, 1e-1, 1e-2]
+# decay: 0.01
\ No newline at end of file
diff --git a/deeplens/geolens.py b/deeplens/geolens.py
index 1822c40..5c4d741 100644
--- a/deeplens/geolens.py
+++ b/deeplens/geolens.py
@@ -28,6 +28,8 @@
 from torchvision.utils import make_grid, save_image
 from tqdm import tqdm
 from transformers import get_cosine_schedule_with_warmup
+import LstsqConvert
+LstsqConvert.register_op()
 
 from .lens import Lens
 from .optics.basics import (
@@ -47,6 +49,7 @@
 from .optics.monte_carlo import forward_integral
 from .optics.ray import Ray
 from .optics.surfaces import (
+    Anamorphic,
     Aperture,
     Aspheric,
     Cubic,
@@ -1061,7 +1064,7 @@ def pupil_field(self, point, wvln=DEFAULT_WAVE, spp=SPP_COHERENT):
         assert spp >= 1000000, (
             "Coherent ray tracing spp is too small, will cause inaccurate simulation."
         )
-        assert torch.get_default_dtype() == torch.float64, (
+        assert torch.get_default_dtype() == torch.float64 or torch.backends.mps.is_available(), (
             "Please set the default dtype to float64 for accurate phase calculation."
         )
 
@@ -1834,6 +1837,39 @@ def calc_entrance_pupil(self, entrance=True, shrink_pupil=False):
             avg_pupilx *= 0.5
         return avg_pupilz, avg_pupilx
 
+    @staticmethod
+    def filter_parallel_lines(Di, Dj, threshold=1e-3):
+        """
+        Identify line pairs that are nearly parallel by computing the cosine of the angle
+        between their directions. If |cos(theta)| ~ 1, they're parallel or anti-parallel.
+
+        Args:
+            Di (torch.Tensor): Direction vectors for line i, shape [M, 2].
+        Dj (torch.Tensor): Direction vectors for line j, shape [M, 2].
+        threshold (float): If 1 - |cos(theta)| < threshold, we consider lines nearly parallel.
+
+        Returns:
+            torch.BoolTensor: A mask of shape [M], where True means "not parallel."
+        """
+        # Dot product for each pair
+        dot_ij = (Di * Dj).sum(dim=-1)  # shape [M]
+        # Norms for each pair
+        norm_i = Di.norm(dim=-1)
+        norm_j = Dj.norm(dim=-1)
+
+        # Avoid division by zero
+        denom = (norm_i * norm_j).clamp_min(1e-12)
+
+        # cos(theta) for each pair
+        cos_theta = dot_ij / denom  # shape [M]
+
+        # Lines are "near parallel" if |cos(theta)| is extremely close to 1
+        # i.e., if 1 - |cos(theta)| < threshold
+        parallel_mask = (1.0 - cos_theta.abs()) < threshold
+
+        # We return the inverse: True means "keep" (not parallel)
+        return ~parallel_mask
+
     @staticmethod
     def compute_intersection_points_2d(origins, directions):
         """Compute the intersection points of 2D lines.
@@ -1846,6 +1882,9 @@ def compute_intersection_points_2d(origins, directions):
             torch.Tensor: Intersection points. Shape: [N*(N-1)/2, 2]
         """
         N = origins.shape[0]
+        if N < 2:
+            # Not enough lines to intersect
+            return torch.empty(0, 2, device=origins.device)
 
         # Create pairwise combinations of indices
         idx = torch.arange(N)
@@ -1856,6 +1895,20 @@ def compute_intersection_points_2d(origins, directions):
         Di = directions[idx_i]  # Shape: [N*(N-1)/2, 2]
         Dj = directions[idx_j]  # Shape: [N*(N-1)/2, 2]
 
+        # 1) Filter out nearly parallel lines
+        not_parallel_mask = GeoLens.filter_parallel_lines(Di, Dj, threshold=1e-3)
+
+        # Keep only the pairs that aren't parallel
+        Oi = Oi[not_parallel_mask]
+        Oj = Oj[not_parallel_mask]
+        Di = Di[not_parallel_mask]
+        Dj = Dj[not_parallel_mask]
+
+        # If everything got filtered out, return empty
+        if Oi.shape[0] == 0:
+            return torch.empty(0, 2, device=origins.device)
+
+        # 2) Construct the system A x = b
         # Vector from Oi to Oj
         b = Oj - Oi  # Shape: [N*(N-1)/2, 2]
 
@@ -1867,7 +1920,10 @@ def compute_intersection_points_2d(origins, directions):
         x, _ = torch.linalg.lstsq(
             A,
             b.unsqueeze(-1),
+            driver='gelsd'
         )[:2]
+        if x.dim() == 2:  # If only one system was solved, add a batch dimension.
+            x = x.unsqueeze(0)
         x = x.squeeze(-1)  # Shape: [N*(N-1)/2, 2]
         s = x[:, 0]
         t = x[:, 1]
@@ -1881,6 +1937,121 @@ def compute_intersection_points_2d(origins, directions):
 
         return P
 
+    def sample_point_source_with_offset(self, depth, num_rays, num_grid, wvln, offset):
+        """
+        Wraps self.sample_point_source to apply a manual offset to the ray source.
+        
+        Args:
+            depth: The object depth.
+            num_rays: Number of rays to sample.
+            num_grid: Grid resolution for sampling.
+            wvln: Wavelength.
+            offset: A tuple (dx, dy) specifying the offset in object space.
+        
+        Returns:
+            A list or batch of rays with modified source positions.
+        """
+        ray = self.sample_point_source(depth=depth, num_rays=num_rays, num_grid=num_grid, wvln=wvln, importance_sampling=True)
+    
+        # Apply offset to ray's origin (or object coordinates).
+        ray.o[..., :2] += torch.tensor(offset, device=self.device)
+        return ray
+
+    def get_cached_rays(self, depth, wvln, offset, num_grid=15):
+        """
+        Retrieve rays sampled at a given depth and wavelength. If not already cached,
+        call sample_point_source once and store the result.
+        
+        Args:
+            depth: The object depth.
+            wvln: Wavelength.
+            offset: A direction specifying the offset in object space. 'c', 'h' or 'v'
+        
+        Returns:
+            A list or batch of rays with modified source positions.
+        """
+        # Create a key based on depth and wavelength
+        key = (depth, wvln, offset, num_grid)
+        if not hasattr(self, "_cached_rays"):
+            self._cached_rays = {}
+        if key not in self._cached_rays:
+            # Sample rays once and store them
+            if offset == 'c':
+                self._cached_rays[key] = self.sample_point_source(
+                    depth=depth,
+                    num_rays=1024,
+                    num_grid=num_grid,
+                    wvln=wvln
+                )
+            elif offset == 'h':
+                self._cached_rays[key] = self.sample_point_source_with_offset(
+                    depth=depth,
+                    num_rays=1024,
+                    num_grid=num_grid,
+                    wvln=wvln,
+                    offset=(1, 0.0)
+                )
+            elif offset == 'v':
+                self._cached_rays[key] = self.sample_point_source_with_offset(
+                    depth=depth,
+                    num_rays=1024,
+                    num_grid=num_grid,
+                    wvln=wvln,
+                    offset=(0.0, 1)
+                )
+        return self._cached_rays[key].clone()
+
+    def compute_effective_horizontal_and_vertical_magnification(self, depth):
+        """
+        Compute the effective horizontal and vertical magnification of the full lens system.
+        This is done by sampling rays from a point source with and without a small
+        horizontal and vertical offset and then comparing the sensor positions.
+    
+        Args:
+            depth (float): The object depth for ray sampling.
+    
+        Returns:
+            torch.Tensor: The computed horizontal and vertical magnification.
+        """
+        delta = 1  # Small horizontal offset in object space.
+        M_x_list = []
+        M_y_list = []
+        for wvln in WAVE_RGB:
+            # Sample a central ray (no offset)
+            ray_central = self.get_cached_rays(
+                depth=depth,
+                wvln=wvln,
+                offset='c'
+            )
+            ray_central, _ = self.trace(ray_central)
+            sensor_central = ray_central.project_to(self.d_sensor)
+
+            # Sample a ray with a small horizontal offset (delta, 0)
+            ray_offset_h = self.get_cached_rays(
+                depth=depth,
+                wvln=wvln,
+                offset='h'
+            )
+            ray_offset_h, _ = self.trace(ray_offset_h)
+            sensor_offset_h = ray_offset_h.project_to(self.d_sensor)
+
+            # Sample a ray with a small vertical offset (0, delta)
+            ray_offset_v = self.get_cached_rays(
+                depth=depth,
+                wvln=wvln,
+                offset='v'
+            )
+            ray_offset_v, _ = self.trace(ray_offset_v)
+            sensor_offset_v = ray_offset_v.project_to(self.d_sensor)
+
+            # Compute horizontal magnification: (difference in sensor x) / delta.
+            M_x = torch.abs((sensor_offset_h[..., 0] - sensor_central[..., 0]) / delta)
+            # Compute vertical magnification: (difference in sensor y) / delta.
+            M_y = torch.abs((sensor_offset_v[..., 1] - sensor_central[..., 1]) / delta)
+            M_x_list.append(M_x)
+            M_y_list.append(M_y)
+        return torch.mean(torch.stack(M_x_list)), torch.mean(torch.stack(M_y_list))
+
     # ====================================================================================
     # Lens operation
     # ====================================================================================
@@ -2773,12 +2944,12 @@ def init_constraints(self):
         else:
             self.is_cellphone = False
 
-            self.dist_min = 0.1
-            self.dist_max = float("inf")
-            self.thickness_min = 0.3
-            self.thickness_max = float("inf")
-            self.flange_min = 0.5
-            self.flange_max = float("inf")
+            self.dist_min = 2.5
+            self.dist_max = 25.0 # float("inf")
+            self.thickness_min = 3.0
+            self.thickness_max = 40.0 # float("inf")
+            self.flange_min = 29.6
+            self.flange_max = 29.8 # float("inf")
 
             self.sag_max = 8.0
             self.grad_max = 1.0
@@ -2812,9 +2983,66 @@ def loss_reg(self, w_focus=None):
                 + 1.0 * loss_surf
                 + 0.05 * loss_angle
             )
+            print(f"Losses in loss_reg: {loss_focus}, {loss_intersec}, {loss_surf}, {loss_angle}")
 
         return loss_reg
 
+    def loss_anamorphic(self, depth, target_local_squeeze, target_global_squeeze):
+        """
+        Compute a combined loss for anamorphic surfaces that includes:
+          - Local penalties on sag, 1st-, and 2nd-order derivatives.
+          - A local squeeze penalty for each anamorphic surface.
+          - A global squeeze penalty computed from ray-traced effective magnification.
+    
+        Args:
+            depth (float): The object depth from which rays are sampled.
+            target_local_squeeze (float): Target local squeeze value.
+            target_global_squeeze (float): Target global squeeze value.
+    
+        Returns:
+            torch.Tensor: Scalar loss value.
+        """
+        total_loss = torch.tensor([0.0], device=self.device)
+
+        # Local loss terms for each anamorphic surface.
+        for idx in self.find_diff_surf():
+            surface = self.surfaces[idx]
+            if not isinstance(surface, Anamorphic):
+                continue
+            # Sample 20 points along the x-axis (from 0 to the aperture radius)
+            x_vals = torch.linspace(0.0, 1.0, 20).to(self.device) * surface.r
+            y_vals = torch.zeros_like(x_vals)
+
+            # --- Sag Penalty ---
+            sag_vals = surface.sag(x_vals, y_vals)
+            # Only penalize if the sag variation exceeds sag_max.
+            sag_penalty = torch.nn.functional.relu((sag_vals.max() - sag_vals.min()) - self.sag_max)
+            total_loss += sag_penalty
+
+            # --- First-order Derivative Penalty ---
+            grad_vals = surface.dfdxyz(x_vals, y_vals)[0]
+            grad_penalty = 10.0 * torch.nn.functional.relu(grad_vals.abs().max() - self.grad_max)
+            total_loss += grad_penalty
+
+            # --- Second-order Derivative Penalty ---
+            grad2_vals = surface.d2fdxyz2(x_vals, y_vals)[0]
+            grad2_penalty = 10.0 * torch.nn.functional.relu(grad2_vals.abs().max() - self.grad2_max)
+            total_loss += grad2_penalty
+            print(f'Total loss: {total_loss}')
+
+        # Global Squeeze Constraint via ray tracing.
+        horiz_mag, vert_mag = self.compute_effective_horizontal_and_vertical_magnification(depth)
+        print(f"Horizontal Magnification: {horiz_mag}")
+        print(f"Vertical Magnification: {vert_mag}")
+        global_squeeze = horiz_mag / (vert_mag + 1e-8)
+        print(f"Global Squeeze: {global_squeeze}")
+        global_squeeze_weight = 1.0 # keeping this big for now - this is a strong constraint
+        global_squeeze_penalty = global_squeeze_weight * (global_squeeze - target_global_squeeze)**2
+        print(f"Global Squeeze Penalty: {global_squeeze_penalty}")
+        total_loss += global_squeeze_penalty
+
+        return total_loss.item()
+
     def loss_infocus(self, bound=0.005):
         """Sample parallel rays and compute RMS loss on the sensor plane, minimize focus loss.
 
@@ -3061,6 +3289,9 @@ def get_optimizer_params(
             elif isinstance(surf, Spheric):
                 params += surf.get_optimizer_params(lr=lr[:2], optim_mat=optim_mat)
 
+            elif isinstance(surf, Anamorphic):
+                params += surf.get_optimizer_params(lr=lr[:2], optim_mat=optim_mat)
+
             else:
                 raise Exception(
                     f"Surface type {surf.__class__.__name__} is not supported for optimization yet."
@@ -3094,6 +3325,7 @@ def optimize(
         match_mat=False,
         shape_control=True,
         importance_sampling=False,
+        anamorphic=False,
         result_dir="./results",
     ):
         """Optimize the lens by minimizing rms errors.
@@ -3107,7 +3339,7 @@ def optimize(
         """
         # Preparation
         depth = DEPTH
-        num_grid = 31
+        num_grid = 15
         spp = 1024
 
         sample_rays_per_iter = 5 * test_per_iter if centroid else test_per_iter
@@ -3175,6 +3407,7 @@ def optimize(
 
             # ===> Optimize lens by minimizing RMS
             loss_rms = []
+            loss_anamorphic_list = []
             for j, wv in enumerate(WAVE_RGB):
                 # Ray tracing
                 ray = rays_backup[j].clone()
@@ -3202,13 +3435,45 @@ def optimize(
                     * weight_mask
                 )
                 loss_rms.append(l_rms)
+                
+                if anamorphic:
+                    # ---- Anamorphic (Squeeze) Loss Computation ----
+                    # Horizontal offset ray.
+                    ray_offset_h = self.get_cached_rays(depth=depth, wvln=wv, offset='h', num_grid=num_grid)
+                    ray_offset_h = self.trace2sensor(ray_offset_h)
+                    xy_offset_h = ray_offset_h.o[..., :2]
+                    ra_xy_offset_h = ray_offset_h.ra.clone().detach()
+                    xy_offset_h_norm = (xy_offset_h - center_p) * ra_xy_offset_h.unsqueeze(-1)
+
+                    # Vertical offset ray.
+                    ray_offset_v = self.get_cached_rays(depth=depth, wvln=wv, offset='v', num_grid=num_grid)
+                    ray_offset_v = self.trace2sensor(ray_offset_v)
+                    xy_offset_v = ray_offset_v.o[..., :2]
+                    ra_xy_offset_v = ray_offset_v.ra.clone().detach()
+                    xy_offset_v_norm = (xy_offset_v - center_p) * ra_xy_offset_v.unsqueeze(-1)
+                    
+                    delta = 1  # object-space offset (ensure consistent units)
+
+                    # Compute effective magnifications.
+                    M_x = torch.mean(torch.abs(xy_offset_h_norm[..., 0])) / delta
+                    M_y = torch.mean(torch.abs(xy_offset_v_norm[..., 1])) / delta
+                    print(f'Magnifications: {M_x}, {M_y}')
+                    ratio = M_x / (M_y + EPSILON)
+                    # Hardcoded for now
+                    # TODO: Make this into a parameter
+                    target_squeeze = 1.5
+                    loss_anamorphic_list.append((ratio - target_squeeze) ** 2)
 
             loss_rms = sum(loss_rms) / len(loss_rms)
+            loss_anamorphic = sum(loss_anamorphic_list) / len(loss_anamorphic_list)
 
             # Total loss
             loss_reg = self.loss_reg()
             w_reg = 0.1
+            w_anamorphic = 0.5
             L_total = loss_rms + w_reg * loss_reg
+            if anamorphic:
+                L_total += w_anamorphic * loss_anamorphic
 
             # Back-propagation
             optimizer.zero_grad()
@@ -3266,6 +3531,9 @@ def read_lens_json(self, filename="./test.json"):
 
                 elif surf_dict["type"] == "ThinLens":
                     s = ThinLens.init_from_dict(surf_dict)
+                
+                elif surf_dict["type"] == "Anamorphic":
+                    s = Anamorphic.init_from_dict(surf_dict)
 
                 else:
                     raise Exception(
@@ -3558,5 +3826,9 @@ def create_surface(surface_type, d_total, aper_r, imgh, mat):
         ai = np.random.randn(7).astype(np.float32) * 1e-30
         k = np.random.randn() * 0.001
         return Aspheric(r=r, d=d_total, c=c, ai=ai, k=k, mat2=mat)
+    elif surface_type == "Plane":
+        return Plane(r=r, d=d_total, mat2=mat)
+    elif surface_type == "Anamorphic":
+        return Anamorphic(r=r, d=d_total, c_x=c, mat2=mat)
     else:
         raise Exception("Surface type not supported yet.")
diff --git a/deeplens/hybridlens.py b/deeplens/hybridlens.py
index 5e971cb..38553fa 100644
--- a/deeplens/hybridlens.py
+++ b/deeplens/hybridlens.py
@@ -45,11 +45,13 @@ class HybridLens(Lens):
     """
     def __init__(self, lens_path):
         super().__init__(lens_path)
-        self.double()
+        if not torch.backends.mps.is_available():
+            self.double()
 
     def double(self):
-        self.geolens.double()
-        self.doe.double()
+        if not torch.backends.mps.is_available():
+            self.geolens.double()
+            self.doe.double()
 
     def read_lens_json(self, lens_path):
         """Read the lens from .json file."""
@@ -233,7 +235,7 @@ def doe_field(self, point, wvln=DEFAULT_WAVE, spp=SPP_COHERENT):
             "which may lead to inaccurate simulation."
         )
         assert (
-            torch.get_default_dtype() == torch.float64
+            torch.get_default_dtype() == torch.float64 or torch.backends.mps.is_available()
         ), "Default dtype must be set to float64 for accurate phase tracing."
 
         geolens, doe = self.geolens, self.doe
@@ -299,7 +301,7 @@ def psf(
             psf_out (torch.Tensor): PSF patch. Normalized to sum to 1. Shape [ks, ks]
         """
         # Check double precision
-        if not torch.get_default_dtype() == torch.float64:
+        if not torch.get_default_dtype() == torch.float64 and not torch.backends.mps.is_available():
             raise ValueError(
                 "Please call HybridLens.double() to set the default dtype to float64 for accurate phase tracing."
             )
diff --git a/deeplens/optics/basics.py b/deeplens/optics/basics.py
index ffb4577..5d17673 100644
--- a/deeplens/optics/basics.py
+++ b/deeplens/optics/basics.py
@@ -11,6 +11,10 @@ def init_device():
         device = torch.device("cuda")
         device_name = torch.cuda.get_device_name(0)
         print(f"Using CUDA: {device_name}")
+    elif torch.backends.mps.is_available():
+        device = torch.device("mps")
+        device_name = "mps"
+        print("Using MPS")
     else:
         device = torch.device("cpu")
         device_name = "CPU"
@@ -264,8 +268,11 @@ def double(self):
 
         torch.set_default_dtype(torch.float64)
         """
+        if torch.backends.mps.is_available():
+            torch.set_default_dtype(torch.float32)
+            return
         assert (
-            torch.get_default_dtype() == torch.float64
+            torch.get_default_dtype() == torch.float64 and not torch.backends.mps.is_available()
         ), "Default dtype should be float64."
 
         for key, val in vars(self).items():
diff --git a/deeplens/optics/ray.py b/deeplens/optics/ray.py
index 3a948da..1e47d29 100644
--- a/deeplens/optics/ray.py
+++ b/deeplens/optics/ray.py
@@ -64,7 +64,7 @@ def propagate_to(self, z, n=1):
         self.o = new_o
 
         if self.coherent:
-            if t.min() > 100 and torch.get_default_dtype() == torch.float32:
+            if t.min() > 100 and torch.get_default_dtype() == torch.float32 and not torch.backends.mps.is_available():
                 raise Warning(
                     "Should use float64 in coherent ray tracing for precision."
                 )
diff --git a/deeplens/optics/surfaces.py b/deeplens/optics/surfaces.py
index e629df0..f2194a3 100644
--- a/deeplens/optics/surfaces.py
+++ b/deeplens/optics/surfaces.py
@@ -2047,3 +2047,281 @@ def surf_dict(self):
         }
 
         return surf_dict
+        
+class Anamorphic(Surface):
+    """Anamorphic surface with different curvatures in x and y directions.
+    
+    This surface implements an anamorphic adapter with different horizontal and
+    vertical curvatures to create a squeeze effect. Typically used in 
+    cinematography to create widescreen images with standard spherical lenses.
+    """
+    
+    def __init__(self, r, d, c_x, mat2, squeeze_factor=1.5, base_focal_length=50.0, is_cylinder=True, device="cpu"):
+        """Initialize an anamorphic surface.
+        
+        Args:
+            r (float): Surface aperture radius
+            d (float or tensor): Distance from origin to surface
+            cx (float): Curvature in x-direction (horizontal)
+            mat2 (str or Material): Material after the surface
+            squeeze_factor (float): Horizontal squeeze factor (default: 1.5x)
+            base_focal_length (float): Base lens focal length in mm (default: 50mm)
+            device (str): Computation device
+        """
+        super(Anamorphic, self).__init__(r, d, mat2, is_square=False, device=device)
+        
+        # Set curvatures for x and y directions
+        self.c_x = torch.tensor(c_x, device=device)
+        self.c_x_perturb = 0.0
+        self.c_y_perturb = 0.0
+        self.is_cylinder = is_cylinder
+        
+        # Store anamorphic parameters
+        self.squeeze_factor = squeeze_factor
+        self.base_focal_length = base_focal_length
+        
+        # Calculate curvatures based on squeeze factor if not provided
+        if c_x is None:
+            # For 1.5x squeeze on 50mm lens:
+            # Vertical focal length = 50mm
+            # Horizontal focal length = 75mm (50mm × 1.5)
+            self.c_x = torch.tensor(1.0 / (base_focal_length * squeeze_factor), device=device)  # Horizontal curvature
+
+        self.to(device)
+    
+    @property
+    def c_y(self):
+        return torch.tensor(1e-9, device=self.device) if self.is_cylinder else self.c_x / self.squeeze_factor
+
+    @c_y.setter
+    def c_y(self, value):
+        if self.is_cylinder:
+            return
+        self.c_x = value * self.squeeze_factor
+
+    @classmethod
+    def init_from_dict(cls, surf_dict):
+        """Initialize surface from a dictionary.
+        
+        Args:
+            surf_dict (dict): Surface parameters dictionary
+            
+        Returns:
+            Anamorphic: Initialized anamorphic surface
+        """
+        if "roc_x" in surf_dict:
+            c_x = 1 / surf_dict["roc_x"]
+        else:
+            c_x = 1 / surf_dict["c_x"]
+        
+        return cls(r=surf_dict["r"], d=surf_dict["d"], c_x=c_x, mat2=surf_dict["mat2"], 
+                  squeeze_factor=surf_dict["squeeze"], base_focal_length=surf_dict["base_focal"])
+        
+    def _sag(self, x, y):
+        """
+        Compute the surface sag for an anamorphic (toroidal) lens.
+        
+        The sag is calculated as the sum of two independent contributions along x and y:
+          sag(x,y) = sag_x(x) + sag_y(y),
+        where:
+          sag_x(x) = (c_x * x^2) / (1 + sqrt(1 - c_x^2 * x^2))
+          sag_y(y) = (c_y * y^2) / (1 + sqrt(1 - c_y^2 * y^2))
+        
+        Here, self.c_x and self.c_y are the nominal curvatures along x and y, and
+        self.c_x_perturb and self.c_y_perturb are any small perturbations.
+        An EPSILON term is added to ensure numerical stability.
+        """
+        c_x = self.c_x + self.c_x_perturb
+        c_y = self.c_y + self.c_y_perturb
+
+        # Handle near-zero curvature to avoid instability:
+        if abs(c_x) < 1e-8:
+            sag_x = torch.zeros_like(x)
+        else:
+            sqrt_arg_x = 1 - c_x**2 * x**2
+            if torch.any(sqrt_arg_x < 0):
+                print("Warning: Negative sqrt_arg_x encountered", sqrt_arg_x)
+            sqrt_arg_x = torch.clamp(sqrt_arg_x, min=EPSILON)
+            A = torch.sqrt(sqrt_arg_x)
+            sag_x = c_x * x**2 / (1 + A)
+
+        if abs(c_y) < 1e-8:
+            sag_y = torch.zeros_like(y)
+        else:
+            sqrt_arg_y = 1 - c_y**2 * y**2
+            if torch.any(sqrt_arg_y < 0):
+                print("Warning: Negative sqrt_arg_y encountered", sqrt_arg_y)
+            sqrt_arg_y = torch.clamp(sqrt_arg_y, min=EPSILON)
+            B = torch.sqrt(sqrt_arg_y)
+            sag_y = c_y * y**2 / (1 + B)
+
+        return sag_x + sag_y
+
+    def _dfdxy(self, x, y):
+        """
+        Compute the first-order derivatives (dz/dx and dz/dy) of the anamorphic sag.
+        Includes a check: if |c_x| or |c_y| is below 1e-8, returns zero derivative.
+        """
+        c_x = self.c_x + self.c_x_perturb
+        c_y = self.c_y + self.c_y_perturb
+
+        # For the x-component:
+        if abs(c_x) < 1e-8:
+            dfdx = torch.zeros_like(x)
+        else:
+            sqrt_arg_x = 1 - c_x**2 * x**2
+            if torch.any(sqrt_arg_x < 0):
+                print("Warning: Negative sqrt_arg_x encountered", sqrt_arg_x)
+            sqrt_arg_x = torch.clamp(sqrt_arg_x, min=EPSILON)
+            A = torch.sqrt(sqrt_arg_x)
+            dfdx = 2 * c_x * x / (1 + A) + (c_x**3 * x**3) / (A * (1 + A)**2)
+
+        # For the y-component:
+        if abs(c_y) < 1e-8:
+            dfdy = torch.zeros_like(y)
+        else:
+            sqrt_arg_y = 1 - c_y**2 * y**2
+            if torch.any(sqrt_arg_y < 0):
+                print("Warning: Negative sqrt_arg_y encountered", sqrt_arg_y)
+            sqrt_arg_y = torch.clamp(sqrt_arg_y, min=EPSILON)
+            B = torch.sqrt(sqrt_arg_y)
+            dfdy = 2 * c_y * y / (1 + B) + (c_y**3 * y**3) / (B * (1 + B)**2)
+
+        return dfdx, dfdy
+
+
+    def _d2fdxy(self, x, y):
+        """
+        Compute the second-order derivatives of the anamorphic sag surface.
+        Returns d²f/dx², d²f/dxdy, d²f/dy².
+        If |c_x| or |c_y| is below 1e-8, returns zero for the corresponding second derivative.
+    
+        Uses a finite difference approximation for the second derivative of the function:
+        F(r²) = c * r² / (1 + sqrt(1 - c² * r²))
+        with r² = x² (or y²).
+        """
+        delta = 1e-6  # Small step for finite difference
+
+        c_x = self.c_x + self.c_x_perturb
+        c_y = self.c_y + self.c_y_perturb
+
+        # For the x-component:
+        if abs(c_x) < 1e-8:
+            d2f_dx2 = torch.zeros_like(x)
+        else:
+            r2_x = x**2
+            S_x = torch.sqrt(torch.clamp(1 - c_x**2 * r2_x, min=EPSILON))
+            Fprime_x = (c_x * (1 + S_x) + c_x**3 * r2_x / (2 * S_x)) / ((1 + S_x)**2)
+            r2_x_plus = r2_x + delta
+            S_x_plus = torch.sqrt(torch.clamp(1 - c_x**2 * r2_x_plus, min=EPSILON))
+            Fprime_x_plus = (c_x * (1 + S_x_plus) + c_x**3 * r2_x_plus / (2 * S_x_plus)) / ((1 + S_x_plus)**2)
+            Fdouble_x = (Fprime_x_plus - Fprime_x) / delta
+
+            d2f_dx2 = 4 * x**2 * Fdouble_x + 2 * Fprime_x
+
+        # For the y-component:
+        if abs(c_y) < 1e-8:
+            d2f_dy2 = torch.zeros_like(y)
+        else:
+            r2_y = y**2
+            S_y = torch.sqrt(torch.clamp(1 - c_y**2 * r2_y, min=EPSILON))
+            Fprime_y = (c_y * (1 + S_y) + c_y**3 * r2_y / (2 * S_y)) / ((1 + S_y)**2)
+            r2_y_plus = r2_y + delta
+            S_y_plus = torch.sqrt(torch.clamp(1 - c_y**2 * r2_y_plus, min=EPSILON))
+            Fprime_y_plus = (c_y * (1 + S_y_plus) + c_y**3 * r2_y_plus / (2 * S_y_plus)) / ((1 + S_y_plus)**2)
+            Fdouble_y = (Fprime_y_plus - Fprime_y) / delta
+
+            d2f_dy2 = 4 * y**2 * Fdouble_y + 2 * Fprime_y
+
+        # Cross derivative is zero since the surface is separable:
+        d2f_dxdy = torch.zeros_like(x)
+
+        return d2f_dx2, d2f_dxdy, d2f_dy2
+
+    def is_within_data_range(self, x, y):
+        """Check if (x, y) is within the valid region of the anamorphic surface.
+
+        For an anamorphic lens the valid region is naturally elliptical. One way to
+        define this is to require that:
+            (cx*x)^2 + (cy*y)^2 < 1
+        """
+        valid_x = (x**2) < 1/(self.c_x**2)
+        valid_y = (y**2) < 1/(self.c_y**2)
+        return valid_x & valid_y
+
+    def max_height(self):
+        """Return the maximum valid y coordinate for the anamorphic surface.
+
+        Since the vertical curvature controls the valid range in y, we use 1/cy.
+        A small margin (e.g., 0.01) is subtracted to avoid boundary issues.
+        """
+        max_y = 1.0 / self.c_y.abs().item() - 0.01
+        return max_y
+
+    def perturb(self, tolerance):
+        """Randomly perturb surface parameters to simulate manufacturing errors.
+
+        In the anamorphic case you might also want to perturb the separate curvatures.
+        """
+        self.r_offset = np.random.randn() * tolerance.get("r", 0.001)
+        self.d_offset = np.random.randn() * tolerance.get("d", 0.001)
+        # Optionally, add curvature perturbations:
+        self.cx_offset = np.random.randn() * tolerance.get("cx", 0.001)
+        # Apply the perturbations if desired:
+        self.c_x = self.c_x + self.cx_offset
+
+    def get_optimizer_params(self, lr=[0.001, 0.001, 0.001], optim_mat=False):
+        """Activate gradient computation for cx, cy, and d and return optimizer parameters.
+
+        Here we include separate parameters for horizontal and vertical curvatures.
+        """
+        self.c_x.requires_grad_(True)
+        self.d.requires_grad_(True)
+
+        params = []
+        params.append({"params": [self.c_x], "lr": lr[0]})
+        params.append({"params": [self.d], "lr": lr[1]})
+
+        if optim_mat and self.mat2.get_name() != "air":
+            params += self.mat2.get_optimizer_params()
+
+        return params
+
+    def surf_dict(self):
+        """Return a dictionary of the anamorphic surface parameters."""
+        roc_x = 1 / self.c_x.item() if self.c_x.item() != 0 else 0.0
+        roc_y = 1 / self.c_y.item() if self.c_y.item() != 0 else 0.0
+        # special: self.r *sometimes* becomes a torch tensor. If it is, unwrap it.
+        r = self.r.item() if isinstance(self.r, torch.Tensor) else self.r
+        surf_dict = {
+            "type": "Anamorphic",
+            "r": round(r, 4),
+            "cx": round(self.c_x.item(), 4),
+            "roc_x": round(roc_x, 4),
+            "cy": round(self.c_y.item(), 4),
+            "roc_y": round(roc_y, 4),
+            "d": round(self.d.item(), 4),
+            "mat2": self.mat2.get_name(),
+        }
+        return surf_dict
+
+    def zmx_str(self, surf_idx, d_next):
+        """Return Zemax surface string for an anamorphic lens."""
+        if self.mat2.get_name() == "air":
+            zmx_str = f"""SURF {surf_idx} 
+    TYPE STANDARD 
+    CURVX {self.c_x.item()} 
+    CURVY {self.c_y.item()} 
+    DISZ {d_next.item()} 
+    DIAM {self.r * 2}
+"""
+        else:
+            zmx_str = f"""SURF {surf_idx} 
+    TYPE STANDARD 
+    CURVX {self.c_x.item()} 
+    CURVY {self.c_y.item()} 
+    DISZ {d_next.item()} 
+    GLAS {self.mat2.get_name().upper()} 0 0 {self.mat2.n} {self.mat2.V}
+    DIAM {self.r * 2}
+"""
+        return zmx_str
diff --git a/lenses/camera/anamorphic_50mm.json b/lenses/camera/anamorphic_50mm.json
new file mode 100644
index 0000000..f4e5079
--- /dev/null
+++ b/lenses/camera/anamorphic_50mm.json
@@ -0,0 +1,221 @@
+{
+    "info": "None",
+    "foclen": 50.0,
+    "fnum": 4.0,
+    "r_sensor": 22.259585847127173,
+    "(d_sensor)": 112.4276,
+    "(sensor_size)": [
+        31.4798,
+        31.4798
+    ],
+    "surfaces": [
+        {
+            "idx": 1,
+            "type": "Anamorphic",
+            "r": 19.4965,
+            "cx": -0.0075,
+            "roc_x": -133.8613,
+            "cy": 0.0,
+            "roc_y": 1000000028.2819,
+            "d": 0.0,
+            "mat2": "1.5785/53.76",
+            "d_next": 6.2061,
+            "squeeze": 1.5,
+            "base_focal": 50.0
+        },
+        {
+            "idx": 2,
+            "type": "Anamorphic",
+            "r": 17.7035,
+            "cx": -0.0164,
+            "roc_x": -61.075,
+            "cy": 0.0,
+            "roc_y": 1000000028.2819,
+            "d": 6.2061,
+            "mat2": "1.6936/56.00",
+            "d_next": 3.0,
+            "squeeze": 1.5,
+            "base_focal": 50.0
+        },
+        {
+            "idx": 3,
+            "type": "Anamorphic",
+            "r": 16.4653,
+            "cx": -0.0077,
+            "roc_x": -129.3221,
+            "cy": 0.0,
+            "roc_y": 1000000028.2819,
+            "d": 9.2061,
+            "mat2": "air",
+            "d_next": 7.7379,
+            "squeeze": 1.5,
+            "base_focal": 50.0
+        },
+        {
+            "idx": 4,
+            "type": "Spheric",
+            "r": 13.6135,
+            "(c)": -0.0042,
+            "roc": -238.2394,
+            "(d)": 16.944,
+            "mat2": "1.5099/63.11",
+            "d_next": 4.5316
+        },
+        {
+            "idx": 5,
+            "type": "Spheric",
+            "r": 11.6499,
+            "(c)": -0.0066,
+            "roc": -152.4719,
+            "(d)": 21.4756,
+            "mat2": "1.5488/59.19",
+            "d_next": 3.0
+        },
+        {
+            "idx": 6,
+            "type": "Spheric",
+            "r": 9.9741,
+            "(c)": -0.001,
+            "roc": -958.4402,
+            "(d)": 24.4756,
+            "mat2": "air",
+            "d_next": 7.7654
+        },
+        {
+            "idx": 7,
+            "type": "Spheric",
+            "r": 7.2474,
+            "(c)": -0.0009,
+            "roc": -1114.5847,
+            "(d)": 32.241,
+            "mat2": "1.6104/44.70",
+            "d_next": 3.0998
+        },
+        {
+            "idx": 8,
+            "type": "Spheric",
+            "r": 6.233,
+            "(c)": -0.0063,
+            "roc": -159.9095,
+            "(d)": 35.3408,
+            "mat2": "air",
+            "d_next": 4.9619
+        },
+        {
+            "idx": 9,
+            "type": "Aperture",
+            "r": 6.25,
+            "(d)": 40.3026,
+            "mat2": "air",
+            "is_square": false,
+            "diffraction": false,
+            "d_next": 4.7149
+        },
+        {
+            "idx": 10,
+            "type": "Spheric",
+            "r": 5.4872,
+            "(c)": -0.0023,
+            "roc": -436.249,
+            "(d)": 45.0175,
+            "mat2": "1.8457/24.85",
+            "d_next": 3.1387
+        },
+        {
+            "idx": 11,
+            "type": "Spheric",
+            "r": 6.6154,
+            "(c)": -0.0115,
+            "roc": -87.0132,
+            "(d)": 48.1561,
+            "mat2": "air",
+            "d_next": 6.3158
+        },
+        {
+            "idx": 12,
+            "type": "Spheric",
+            "r": 8.3801,
+            "(c)": -0.0117,
+            "roc": -85.7145,
+            "(d)": 54.472,
+            "mat2": "1.5300/56.09",
+            "d_next": 3.1915
+        },
+        {
+            "idx": 13,
+            "type": "Spheric",
+            "r": 9.3873,
+            "(c)": -0.0188,
+            "roc": -53.1772,
+            "(d)": 57.6634,
+            "mat2": "1.5169/65.29",
+            "d_next": 3.4949
+        },
+        {
+            "idx": 14,
+            "type": "Spheric",
+            "r": 10.2942,
+            "(c)": -0.0296,
+            "roc": -33.7922,
+            "(d)": 61.1584,
+            "mat2": "air",
+            "d_next": 8.0985
+        },
+        {
+            "idx": 15,
+            "type": "Anamorphic",
+            "r": 12.5272,
+            "cx": -0.0062,
+            "roc_x": -160.1566,
+            "cy": 0.0,
+            "roc_y": 1000000028.2819,
+            "d": 69.2568,
+            "mat2": "1.6935/55.78",
+            "d_next": 2.9999,
+            "squeeze": 1.5,
+            "base_focal": 50.0
+        },
+        {
+            "idx": 16,
+            "type": "Anamorphic",
+            "r": 13.1711,
+            "cx": -0.0031,
+            "roc_x": -318.8873,
+            "cy": 0.0,
+            "roc_y": 1000000028.2819,
+            "d": 72.2567,
+            "mat2": "air",
+            "d_next": 7.5709,
+            "squeeze": 1.5,
+            "base_focal": 50.0
+        },
+        {
+            "idx": 17,
+            "type": "Anamorphic",
+            "r": 15.5294,
+            "cx": -0.0026,
+            "roc_x": -385.9604,
+            "cy": 0.0,
+            "roc_y": 1000000028.2819,
+            "d": 79.8276,
+            "mat2": "1.6936/55.97",
+            "d_next": 3.0,
+            "squeeze": 1.5,
+            "base_focal": 50.0
+        },
+        {
+            "idx": 18,
+            "type": "Anamorphic",
+            "r": 16.3245,
+            "cx": -0.0003,
+            "roc_x": -3447.7029,
+            "cy": 0.0,
+            "roc_y": 1000000028.2819,
+            "d": 82.8276,
+            "mat2": "air",
+            "d_next": 29.6,
+            "squeeze": 1.5,
+            "base_focal": 50.0
+        }
+    ]
+}
\ No newline at end of file
diff --git a/lenses/camera/anamorphic_custom_f2.0.json b/lenses/camera/anamorphic_custom_f2.0.json
new file mode 100644
index 0000000..dd7bc29
--- /dev/null
+++ b/lenses/camera/anamorphic_custom_f2.0.json
@@ -0,0 +1,57 @@
+{
+    "info": "None",
+    "foclen": 50.0,
+    "fnum": 1.8,
+    "r_sensor": 21.6,
+    "d_sensor": 109.25,
+    "sensor_size": [
+        30.547012947258853,
+        30.547012947258853
+    ],
+    "surfaces": [
+        {
+            "type": "Plane",
+            "l": 37.4766594029,
+            "r": 26.5,
+            "d": 0.0,
+            "is_square": true,
+            "mat2": "n-sf10",
+            "d_next":  6.73
+        },
+        {
+            "type": "Anamorphic",
+            "r": 26.5,
+            "c_x": 0.02666666667,
+            "roc_x": 37.5,
+            "c_y": 0.000000001,
+            "roc_y": 1000000000,
+            "squeeze": 1.5,
+            "base_focal": 50.0,
+            "d": 6.73,
+            "mat2": "air",
+            "d_next": 17.77
+        },
+        {
+            "type": "Anamorphic",
+            "r": 17.5,
+            "c_x": 0.02,
+            "roc_x": 50.0,
+            "c_y": 0.000000001,
+            "roc_y": 1000000000,
+            "squeeze": 1.5,
+            "base_focal": 50.0,
+            "d": 24.5,
+            "mat2": "n-sf10",
+            "d_next": 4.0
+        },
+        {
+            "type": "Plane",
+            "l": 24.7487373415,
+            "r": 17.5,
+            "d": 28.5,
+            "is_square": true,
+            "mat2": "air",
+            "d_next": 80.75
+        }
+    ]
+}
diff --git a/lenses/camera/iter17500.json b/lenses/camera/iter17500.json
new file mode 100644
index 0000000..ae67ecd
--- /dev/null
+++ b/lenses/camera/iter17500.json
@@ -0,0 +1,63 @@
+{
+    "info": "None",
+    "foclen": 37.5,
+    "fnum": 2.0,
+    "r_sensor": 40.5,
+    "(d_sensor)": 109.0547,
+    "(sensor_size)": [
+        30.4056,
+        30.4056
+    ],
+    "surfaces": [
+        {
+            "idx": 1,
+            "type": "Plane",
+            "(l)": 34.64823227814083,
+            "r": 24.5,
+            "(d)": 0.0,
+            "is_square": true,
+            "mat2": "n-sf10",
+            "d_next": 3.5667
+        },
+        {
+            "idx": 2,
+            "type": "Spheric",
+            "r": 24.5,
+            "(c)": 0.008695652174,
+            "roc": 115.0,
+            "(d)": 3.5667,
+            "mat2": "air",
+            "d_next": 18.8333
+        },
+        {
+            "idx": 3,
+            "type": "Aperture",
+            "r": 17.7,
+            "(d)": 22.4,
+            "mat2": "air",
+            "is_square": false,
+            "diffraction": false,
+            "d_next": 0.1
+        },
+        {
+            "idx": 4,
+            "type": "Spheric",
+            "r": 17.5,
+            "(c)": 0.0164,
+            "roc": 61.1577,
+            "(d)": 22.5,
+            "mat2": "n-sf10",
+            "d_next": 3.165
+        },
+        {
+            "idx": 5,
+            "type": "Plane",
+            "(l)": 24.7487373415,
+            "r": 17.5,
+            "(d)": 25.665,
+            "is_square": true,
+            "mat2": "air",
+            "d_next": 83.3897
+        }
+    ]
+}
\ No newline at end of file