(* 交互演示：范数可视仪 —— RMS vs PV *)
(* 理解不同的数学"尺子"如何评价光学误差 *)


Manipulate[
 Module[{xRange, idealProfile, errorProfile, actualProfile, norm1, norm2, normInf, errorText},
  
  xRange = Range[-1, 1, 0.01];
  
  (* 1. 定义完美轮廓 (理想的球面) *)
  idealProfile = 1 - xRange^2;
  
  (* 2. 制造误差 *)
  errorProfile = Switch[errorType,
    "整体偏移 (Defocus)", 
      ConstantArray[magnitude, Length[xRange]], (* 整个镜片度数做错了 *)
      
    "局部尖峰 (Spike)", 
      (* 镜片中间有个凸起小点，模拟灰尘或加工瑕疵 *)
      magnitude * Exp[-100 * xRange^2], 
      
    "随机噪点 (Roughness)", 
      (* 表面粗糙度 *)
      SeedRandom[1];
      magnitude * RandomReal[{-1, 1}, Length[xRange]]
  ];
  
  actualProfile = idealProfile + errorProfile;
  
  (* 3. 计算三种范数 (三把尺子) *)
  (* L1: 平均绝对误差 *)
  norm1 = Norm[errorProfile, 1] / Length[xRange]; 
  (* L2: 均方根误差 (RMS) *)
  norm2 = Norm[errorProfile, 2] / Sqrt[Length[xRange]]; 
  (* L_inf: 最大误差 (Peak Error) *)
  normInf = Norm[errorProfile, Infinity];
  
  (* 4. 可视化 *)
  Column[{
    (* 上图：镜片轮廓对比 *)
    Show[
      Plot[1 - t^2, {t, -1, 1}, PlotStyle -> {Thick, Blue}, PlotLegends -> {"理想设计 (Target)"}],
      ListLinePlot[Transpose[{xRange, actualProfile}], PlotStyle -> {Orange}, PlotLegends -> {"实际加工 (Actual)"}],
      PlotRange -> All, ImageSize -> 450, AspectRatio -> 1/3,
      PlotLabel -> Style["镜片截面图", Bold, 14],
      Filling -> {2 -> {1}}, FillingStyle -> Opacity[0.2, Red]
    ],
    
    (* 下图：误差分析仪表盘 *)
    Grid[{
      {Style["数学尺子 (范数)", Bold], Style["数值", Bold], Style["临床意义", Bold]},
      
      {"L1 范数 (Mean Abs)", 
       NumberForm[norm1, {1, 3}], 
       "平均厚度偏差"},
      
      {Style["L2 范数 (RMS)", Blue, Bold], 
       Style[NumberForm[norm2, {1, 3}], Blue, Bold], 
       Style["整体像质 (Spot Size)", Blue]},
       
      {Style["L∞ 范数 (Max/PV)", Red, Bold], 
       Style[NumberForm[normInf, {1, 3}], Red, Bold], 
       Style["最坏点 (最大像差)", Red]}
    }, Frame -> All, Background -> {{LightGray, None}, None}]
  }]
 ],
 
 (* 控制面板 *)
 Style["误差模拟器", Bold, 12],
 {{errorType, "局部尖峰 (Spike)", "误差类型"}, {"整体偏移 (Defocus)", "局部尖峰 (Spike)", "随机噪点 (Roughness)"}},
 {{magnitude, 0.2, "误差幅度"}, 0, 0.5, Appearance -> "Labeled"}
]

(* 代码处方 08 (Debug版)：镜片良品率的蒙特卡洛模拟 *)

(* 0. 清除之前的定义，防止符号冲突 *)
Clear["Global`*"]; 

(* 1. 定义设计目标 *)
targetPower = -3.00; (* 目标度数 *)
tolerance   = 0.12;  (* 允许误差 *)

(* 2. 定义加工误差模型 *)
(* 【关键修正】：变量名去掉下划线，改用驼峰命名 *)

(* 设备 A: 正态分布 (Normal Distribution) *)
machineASigma = 0.08; 

(* 设备 B: 学生 t-分布 (Student's t-Distribution) *)
(* 修正说明：Wolfram 中下划线是保留字符，machineB_Scale 写法是非法的 *)
machineBScale = 0.03; 
machineBDf    = 3;    

(* 3. 执行蒙特卡洛模拟 *)
numberOfLenses = 2000;
Print["正在模拟生产 ", numberOfLenses, " 片镜片..."];

(* 生产批次 A：直接生成 *)
batchA = RandomVariate[NormalDistribution[targetPower, machineASigma], numberOfLenses];

(* 生产批次 B：手动构建位置-尺度 T 分布 *)
(* 公式：数据 = 目标值 + 尺度 * 标准T分布随机数 *)
batchB = targetPower + machineBScale * RandomVariate[StudentTDistribution[machineBDf], numberOfLenses];

(* 4. 统计良品率 (向量化计算) *)
(* 使用 Count 直接统计符合条件的个数，效率极高 *)
countPassA = Count[batchA, x_ /; Abs[x - targetPower] <= tolerance];
countPassB = Count[batchB, x_ /; Abs[x - targetPower] <= tolerance];

passRateA = N[countPassA / numberOfLenses * 100];
passRateB = N[countPassB / numberOfLenses * 100];

(* 5. 可视化对比报告 *)
Print["\n================================="];
Print[Style["公差分析报告 (Debug 完成版)", 18, Bold, FontFamily -> "Arial"]];
Print["================================="];
Print["目标度数: ", targetPower, " D | 允许误差: +/- ", tolerance, " D"];
Print["---------------------------------"];
Print["设备 A (低精度稳定型) 良品率: ", Style[Row[{passRateA, "%"}], Bold, Blue]];
Print["设备 B (高精度长尾型) 良品率: ", Style[Row[{passRateB, "%"}], Bold, Orange]];

(* 智能结论 *)
If[passRateA > passRateB,
    Print[Style["\n结论：设备 A 胜出。设备 B 虽然精度高(Scale小)，但由于 T分布的长尾特性(df=3)，导致出现了更多严重超标的废品。", Gray]],
    Print[Style["\n结论：设备 B 胜出。", Gray]]
];

(* 绘制直方图对比 *)
PairedHistogram[batchA, batchB, 
    50, 
    "PDF", 
    ChartStyle -> {Blue, Orange},
    PlotLabel -> Style["度数分布对比：A (正态) vs B (长尾)", 14, Bold],
    AxesLabel -> {"度数 (D)", "概率密度"},
    PlotRange -> {{targetPower - 0.5, targetPower + 0.5}, All}, (* 调整范围以便看清中心 *)
    
    (* 红色虚线标出公差范围 *)
    GridLines -> {{targetPower - tolerance, targetPower + tolerance}, None},
    GridLinesStyle -> Directive[Red, Thick, Dashed],
    
    ChartLegends -> {"设备 A (Normal)", "设备 B (Student-t)"}
]

(* 交互演示：寻找帕累托前沿 (Pareto Frontier) *)
(* 模拟光学设计中的 Trade-off：想要中心度数足 (Focus)，又想周边像差小 (Flatness) *)


Manipulate[
 Module[{x, shape, focusScore, flatnessScore, weightedObj, optSol, shapeFunc},
  
  (* 1. 定义我们的镜片模型：一个简单的多项式 f(x) = a*x^2 + b*x^4 *)
  (* 我们限制在 x=[-1, 1] 区间 *)
  shape[x_, a_] := a * (1 - x^2)^2; 
  
  (* 2. 定义两个冲突的目标 *)
  (* 目标 A：中心屈光力要足 (maximize height at x=0) *)
  (* 我们用负号，因为 NMinimize 是求极小值 *)
  (* Target: shape(0) 越接近 1.0 越好 *)
  focusScore[a_] := (shape[0, a] - 1.0)^2;
  
  (* 目标 B：表面要平坦 (minimize slope/distortion) *)
  (* 计算 x=0 到 1 的平均斜率平方 *)
  flatnessScore[a_] := Integrate[D[shape[x, a], x]^2, {x, 0, 1}];
  
  (* 3. 标量化 (Scalarization)：将两个目标合成一个 *)
  (* Total Cost = w * Objective1 + (1-w) * Objective2 *)
  weightedObj[a_] := weight * focusScore[a] + (1 - weight) * flatnessScore[a];
  
  (* 4. 求解最优的参数 a *)
  optSol = NMinimize[{weightedObj[a], 0 <= a <= 1.5}, a];
  bestA = a /. optSol[[2]];
  
  (* 5. 计算当前两个具体的分数 *)
  currentFocusError = focusScore[bestA];
  currentFlatnessError = flatnessScore[bestA];
  
  (* 6. 可视化 *)
  Grid[{{
    (* 左图：镜片形状展示 *)
    Plot[shape[x, bestA], {x, -1, 1},
     PlotRange -> {-0.1, 1.2}, ImageSize -> 350,
     PlotStyle -> {Thick, Purple},
     Filling -> Axis, FillingStyle -> Opacity[0.1, Purple],
     PlotLabel -> Style["镜片截面形状", Bold, 14],
     Epilog -> {
       Text["中心高度 (清晰度)", {0, shape[0, bestA]}, {0, -1.5}],
       Text["边缘坡度 (像差)", {0.8, 0.5}, {-1, 0}]
     }],
     
    (* 右图：帕累托前沿图 *)
    Show[
      (* 预先计算好的帕累托曲线背景 (近似) *)
      ParametricPlot[{ (a-1)^2, Integrate[(D[a*(1-t^2)^2,t])^2,{t,0,1}] }, {a, 0, 1.2},
       PlotStyle -> {Dashed, Gray}, AspectRatio -> 1, PlotRange -> All],
       
      (* 当前的工作点 *)
      Graphics[{
        PointSize[0.03], Red, Point[{currentFocusError, currentFlatnessError}],
        Text[Style["当前设计点", Red, Bold], {currentFocusError, currentFlatnessError}, {-1.2, 0}]
      }],
      
      AxesLabel -> {"视力误差 (Focus Error)", "像差水平 (Distortion)"},
      PlotLabel -> Style["帕累托前沿 (Pareto Frontier)", Bold, 14],
      ImageSize -> 350,
      Epilog -> {Inset[Style["越靠近原点越好\n但无法进入虚线内", 10, Gray], {0.5, 0.5}]}
    ]
  }}]
 ],
 
 (* 控制面板 *)
 Style["设计权衡控制", 12, Bold],
 {{weight, 0.5, "偏好权重 (w)"}, 0.01, 0.99, Appearance -> "Labeled", 
  ControlPlacement -> Top},
 Row[{Style["<-- 重视舒适 (Soft)", Darker[Green]], Spacer[150], Style["重视清晰 (Hard) -->", Blue]}]
]

(* 代码处方 09：生成帕累托前沿分析报告 *)


(* 1. 定义模拟的目标函数 (Cost Functions) *)
(* 这里使用抽象函数模拟复杂的光学光线追踪结果 *)
(* f1: 清晰度损失 (越小越好) *)
costClear[x_] := (x - 2)^2; 
(* f2: 舒适度损失 (越小越好) *)
costComfort[x_] := (x + 1)^2 + 2; 


(* 注意：这两个函数最小值位置不同。f1想去 x=2, f2想去 x=-1 *)


(* 2. 扫描权重，生成前沿数据 *)
points = Table[
   Module[{w, sol, xOpt, val1, val2},
    w = weight;
    
    (* 求解加权总分最小化 *)
    (* Constrained optimization *)
    sol = NMinimize[{w * costClear[x] + (1 - w) * costComfort[x], -2 <= x <= 3}, x];
    xOpt = x /. sol[[2]];
    
    (* 记录两个独立目标的得分 *)
    val1 = costClear[xOpt];
    val2 = costComfort[xOpt];
    
    (* 返回数据包: {权重, 视力损失, 舒适损失, 最优参数x} *)
    {w, val1, val2, xOpt}
   ],
   {weight, 0, 1, 0.05} (* 从 0 到 1，步长 0.05 *)
   ];


(* 3. 数据提取与处理 *)
paretoFrontier = points[[All, {2, 3}]]; (* 提取坐标点 {J1, J2} *)


(* 4. 绘制帕累托前沿图 *)
Show[
 ListPlot[paretoFrontier, 
  PlotStyle -> {PointSize[0.02], Darker[Blue]},
  Joined -> True, InterpolationOrder -> 2,
  PlotTheme -> "Scientific",
  FrameLabel -> {"清晰度损失 (Clarity Loss)", "舒适度损失 (Comfort Loss)"},
  PlotLabel -> Style["设计方案的帕累托前沿", 16, Bold],
  GridLines -> Automatic,
  Epilog -> {
    Text["硬性设计\n(Hard Design)", points[[-1, {2, 3}]] + {-0.5, 0.5}],
    Text["软性设计\n(Soft Design)", points[[1, {2, 3}]] + {0.5, -0.5}],
    Text["平衡设计\n(Balanced)", points[[10, {2, 3}]] + {0.5, 0.5}],
    Red, PointSize[0.03], Point[points[[10, {2, 3}]]]
  }
 ],
 Graphics[{Dashed, Gray, Line[{{0,0}, {10,10}}], Text["理想原点(不可达)", {0.5, 0.5}]}]
]


(* 5. 输出决策建议表 *)
Print["--------------------------------------------------"];
Print[Style["设计方案建议表 (Top 3)", 14, Bold]];
Print["--------------------------------------------------"];
Print[Grid[{
   {"方案类型", "偏好权重(w)", "清晰度评分", "舒适度评分", "适用人群"},
   {"软性设计 (Soft)", NumberForm[points[[2, 1]], {2, 2}], NumberForm[10 - points[[2, 2]], {2, 1}], NumberForm[10 - points[[2, 3]], {2, 1}], "初次配戴/敏感者"},
   {"平衡设计 (Balanced)", NumberForm[points[[11, 1]], {2, 2}], NumberForm[10 - points[[11, 2]], {2, 1}], NumberForm[10 - points[[11, 3]], {2, 1}], "大多数用户"},
   {"硬性设计 (Hard)", NumberForm[points[[-2, 1]], {2, 2}], NumberForm[10 - points[[-2, 2]], {2, 1}], NumberForm[10 - points[[-2, 3]], {2, 1}], "长期配戴/静止工作"}
  }, Frame -> All, Background -> {{LightGray, None}, None}]];

(* 交互演示：视觉模拟器 —— 卷积的力量 *)
(* 模拟不同像差形式的 PSF 如何"涂抹"视力表 *)


Manipulate[
 Module[{sourceImage, psfKernel, simulation},
  
  (* 1. 生成源图像：一个标准的 Snellen E 视标 *)
  sourceImage = Rasterize[
    Graphics[{Black, Rectangle[{-2, -2}, {2, 2}], White, 
      Text[Style["E", FontFamily -> "Helvetica", FontSize -> 80, Bold], {0, 0}]},
     PlotRange -> {{-1, 1}, {-1, 1}}],
    ImageSize -> 200, RasterSize -> 128, ColorSpace -> "Grayscale"];
  
  (* 2. 构建点扩散函数 (PSF) / 卷积核 *)
  (* 这是一个二维高斯分布，我们可以拉伸它模拟散光 *)
  psfKernel = Table[
    Exp[-( (x^2)/(2*sigma^2) + (y^2)/(2*(sigma*ratio)^2) )],
    {y, -3, 3, 0.2}, {x, -3, 3, 0.2}
  ];
  
  (* 旋转 PSF 模拟散光轴位 *)
  psfKernel = ImageRotate[Image[psfKernel, "Real"], axis Degree, Background -> 0];
  
  (* 3. 执行卷积 (Convolution) *)
  (* 核心函数：ImageConvolve *)
  (* 这模拟了光学系统对图像的"传递"过程 *)
  simulation = ImageConvolve[sourceImage, psfKernel];
  
  (* 4. 可视化 *)
  Column[{
    Row[{
      Column[{Style["1. 原始物体", Bold], sourceImage}, Alignment -> Center],
      Style[" * ", FontSize -> 30], (* 卷积符号 *)
      Column[{Style["2. 眼睛的 PSF", Bold], ImageResize[psfKernel, 50]}, Alignment -> Center],
      Style[" = ", FontSize -> 30],
      Column[{Style["3. 患者看到的", Bold, Red], simulation}, Alignment -> Center]
    }, Alignment -> Center],
    
    Spacer[20],
    Style[StringTemplate["散光轴位: ``° | 散光量: ``"][NumberForm[axis, {3, 0}], NumberForm[ratio, {2, 1}]], Gray]
  }, Alignment -> Center]
 ],
 
 (* 控制面板 *)
 Style["像差参数控制", 12, Bold],
 {{sigma, 0.5, "模糊程度 (Defocus)"}, 0.1, 1.5},
 {{ratio, 1.0, "散光拉伸 (Cylinder)"}, 1.0, 5.0},
 {{axis, 0, "散光轴位 (Axis)"}, 0, 180, 1}
]

(* 代码处方 10 (修正版)：衍射极限 PSF 模拟器 *)
(* 利用 FFT 计算光瞳形状对成像质量的影响 *)

(* ---------------------------------------------------- *)
(* 1. 定义光瞳函数 (Pupil Function) *)
(* ---------------------------------------------------- *)
n = 512; (* 提高分辨率到 512 以获得更细腻的艾里斑 *)
L = 2.0; (* 物理尺寸范围 [-1, 1] *)

(* 技巧：直接生成坐标矩阵，比 CoordinateBoundsArray 后 Map 效率高百倍 *)
xRange = Subdivide[-1., 1., n - 1];
yRange = Subdivide[-1., 1., n - 1];

(* 使用 Outer 生成网格坐标的模长平方 (r^2 = x^2 + y^2) *)
(* 这是完全向量化的操作，无需 Map *)
rSquaredMatrix = Outer[Plus, xRange^2, yRange^2];

(* 定义圆形光瞳 (Aperture) *)
pupilRadius = 0.2; (* 瞳孔半径 *)
(* UnitStep[t]：当 t >= 0 时为 1，否则为 0。比 If 更快 *)
aperture = UnitStep[pupilRadius^2 - rSquaredMatrix]; 


(* ---------------------------------------------------- *)
(* 2. 计算傅里叶变换 (The Fourier Transform) *)
(* ---------------------------------------------------- *)
(* 物理魔法：光瞳的远场衍射图样 = 光瞳的傅里叶变换 *)

(* FourierParameters -> {1, 1} 是现代光学的常用定义 *)
(* 但默认参数 {0, 1} 对于计算强度分布(模平方)也是完全正确的 *)
wavefront = Fourier[aperture];

(* 修正点：使用 RotateRight 实现标准的 fftshift *)
(* 将零频率分量从数组角落(1,1)移动到数组中心 *)
centeredWavefront = RotateRight[wavefront, {n/2, n/2}];


(* ---------------------------------------------------- *)
(* 3. 获得点扩散函数 (PSF) *)
(* ---------------------------------------------------- *)
(* PSF 是振幅的平方 *)
psf = Abs[centeredWavefront]^2;

(* 为了模拟人眼/相机的对数响应，同时避免中心过曝 *)
(* 加上 1.0 是为了防止 Log[0] *)
logPSF = Log[1.0 + psf];


(* ---------------------------------------------------- *)
(* 4. 可视化对比 (重点改进) *)
(* ---------------------------------------------------- *)

(* 裁剪：为了看清中心的衍射环，我们只取中心 100x100 的区域 *)
(* 直接切片矩阵是最稳健的"放大镜" *)
zoomRadius = 50; 
centerIndex = n/2;
(* Span (;;) 语法：从 center - 50 到 center + 50 *)
zoomedPSF = logPSF[[ centerIndex - zoomRadius ;; centerIndex + zoomRadius, 
                     centerIndex - zoomRadius ;; centerIndex + zoomRadius ]];

GraphicsRow[{
  (* 图1：瞳孔形状 *)
  ArrayPlot[aperture, 
    ColorFunction -> "GrayTones", 
    Frame -> False, 
    PlotLabel -> Style["1. 瞳孔形状 (Aperture)", 14, Bold, FontFamily->"Arial"]],
  
  (* 中间：变换过程 *)
  Graphics[{
     Text[Style["FFT\nSquared", 12, Gray], {0, 0}], 
     Arrow[{{-1, -0.5}, {1, -0.5}}]
  }, ImageSize -> {50, 50}],
  
  (* 图2：衍射斑 (Airy Disk) *)
  ArrayPlot[zoomedPSF, 
    ColorFunction -> "SunsetColors", (* 这种配色能清晰显示次级光环 *)
    Frame -> False,
    PlotLabel -> Style["2. 视网膜成像 (PSF)", 14, Bold, FontFamily->"Arial"],
    Epilog -> {
        (* 加上标注圈出第一个暗环 *)
        {White, Dashed, Thin, Circle[{zoomRadius + 1, zoomRadius + 1}, 6]}
    }
   ]
}, ImageSize -> 700]


(* ---------------------------------------------------- *)
(* 5. 临床与物理思考 *)
(* ---------------------------------------------------- *)
Print["\n================ 代码解读 ================"];
Print["1. 中心亮斑 (Airy Disk)：承载了图像的主要信息。亮斑越小，视力越好（分辨率越高）。"];
Print["2. 衍射环 (Rings)：注意右图中中心红点周围的淡淡光环。"];
Print["   这是光的波动性带来的物理限制。光不再是一个完美的'点'，而是一个'斑'。"];
Print["3. 临床意义：如果瞳孔不仅是圆的，还有多焦晶体的衍射台阶，"];
Print["   这些光环的能量会增强，导致患者在夜间看车灯时出现明显的'光晕' (Halo)。"];