泽尼克多项式

clc;
clear all;
format long
%R0;被测镜半径
%Dia;被测镜口径
Dia=10;
x_Number=401;%x方向采样数
y_Number=x_Number;%y方向采样数
%qiu_o=[00-R0];球心坐标 
%H0=（R0^2-（Dia/2）^2）^0.5;

%% 在被测镜面上生成点阵
%[x_lensy_lens]=meshgrid（-Dia/2:Dia/（x_Number-1）:Dia/2Dia/2:-Dia/（y_Number-1）:-Dia/2）;
%for i=1:1:x_Number;
  %for j=1:1:y_Number;  
%      if （x_lens（ij）^2+y_lens（ij）^2）^0.5>Dia/2
%          x_lens（ij）=NaN;
%          y_lens（ij）=NaN;
%      end
%  end
%end

[x_lensy_lens] = meshgrid（-1:0.01:1-1:0.01:1）;
x=x_lens;
y=y_lens;
WF = zeros（size（x））+nan;
dzdx = zeros（size（x））+nan;
dzdy = zeros（size（x））+nan;
ind = find（x.^2+y.^2 <= 1）;
x1 = x（ind）;
y1 = y（ind）;


Cc（1:36）=1;
%Cc（1）=1;Cc（2）=1;Cc（3）=1;Cc（4）=1;Cc（5）=1;Cc（6）=1;Cc（7）=1;Cc（8）=1;Cc（9）=1;Cc（10）=1;Cc（11）=1;Cc（12）=1;

        z_fitting=Cc（1）*x+...
        Cc（2）*y+...
        Cc（3）*（2*x.^2+2*y.^2-1）+...
        Cc（4）*（x.^2-y.^2）+...
        Cc（5）*2*x.*y+...
        Cc（6）*（3*x.^3+3*x.*y.^2-2*x）+...
        Cc（7）*（3*x.^ 2.*y+3*y.^3-2*y）+...
        Cc（8）*（6*x.^4-6*x.^2+1）+...
        Cc（9）*（x.^3-3*x.*y.^2）+...
        Cc（10）*（3*x.^2.*y-y.^3）+...
        Cc（11）*（4*x.^4-3*x.^2+3*y.^2-4*y.^4）+...
        Cc（12）*（8*x.^3.*y+8*x.*y.^3-6*x.*y）+...
        Cc（13）*（10*x.^5+20*x.^3.*y.^2+10*x.*y.^4-12*x.^3-12*x.*y.^2+3*x）+...
        Cc（14）*（10*x.^4.*y+20*x.^2.*y.^3+10*y.^5-12*x.^2.*y-12*y.^3+3*y）+...
        Cc（15）*（20*x.^6-30*x.^4+12*x.^2-1）+...
        Cc（16）*（x.^4-6*x.^2.*y.^2+y.^4）+...
        Cc（17）*（4*x.^3.*y-4*x.*y.^3）+...
        Cc（18）*（5*x.^5-10*x.^3.*y.^2-4*x.^3-15*x.*y.^4+12*x.*y.^2）+...
        Cc（19）*（-5*y.^5+10*x.^2.*y.^3+4*y.^3+15*y.*x.^4-12*y.*x.^2）+...
        Cc（20）*（15*x.^6+15*x.^4.*y.^2-15*x.^2.*y.^4+6*x.^2-6*y.^2-15*y.^6-20*x.^4+20*y.^4）+...
        Cc（21）*（30*x.^5.*y+60*x.^3.*y.^3+30*x.*y.^5-40*x.^3.*y-40*x.*y.^3+12*x.*y）+...
        Cc（22）*（35*x.^7+105*x.^5.*y.^2+105*x.^3.*y.^4+35*x.*y.^6-60*x.^5-120*x.^3.*y.^2-60*x.*y.^4+30*x.^3+30*x.