用OxyPlot在WPF中演示正演磁异常的变化规律

为了在展示实验成果时动态演示理论球体磁异常随其埋深、磁化倾角的变化规律，我用WPF写了一个小程序来作演示。

MatLab计算磁异常数据

首先是计算理论球体磁异常数据，在Matlab中可以很方便地计算。
为了在展示时能够同时改变磁化倾角（is)、埋深（Deep），计算数据时套了两层循环，外层是is从0到90°，内层是Deep从10到100m，里面计算一条剖面上的Za和Hax磁异常数据。实现代码如下：

function Za=sphere_deeps()
    % 测点分布范围
    dx=5; % X方向测点间距
    nx=81; % X方向测点数
    xmin=-200; % X方向起点
    x=xmin:dx:(xmin+(nx-1)*dx); % X方向范围

    u=4*pi*10^(-7);  %磁导率
    i=pi/2;  %有效磁化倾角is
    a=0; %剖面磁方位角
    T=50000;%地磁场T=50000nT

    % 球体参数
    R1=10; % 球体半径 m
    %D1=30; % 球体埋深 m
    v1=4*pi*R1^3;
    k=0.2; %磁化率
    M1=k*T/u;  %磁化强度 A/m
    m1=M1*v1;   %磁矩

    % 球体 理论磁异常主剖面
    %(x-0),(y-0)
    y=0;
    Za=zeros(46,nx);
    Hax=zeros(46,nx);

    fp=fopen('za_all.out','w');
    for ii=0:90
        i=ii*pi/180;
        for j=0:45
            D1=10+j*2;
            Za0=(u*m1*((2*D1.^2-x.^2-y.^2)*sin(i)-3*D1*x.*cos(i)*cos(a)-3*D1*y.*cos(i)*sin(a)))./(4*pi*(x.^2+y.^2+D1.^2).^(5/2));
            for k=1:nx
                fprintf(fp,'%g ',Za0(k));
            end
            fprintf(fp,'\n');
        end
    end
    fclose(fp);

    fp=fopen('hax_all.out','w');
    for ii=0:90
        i=ii*pi/180;
        for j=0: