为了在展示实验成果时动态演示理论球体磁异常随其埋深、磁化倾角的变化规律，我用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:45
D1=10+j*2;
Hax0=(u*m1*((2*x.^2-y.^2-D1.^2)*cos(i)*cos(a)-3*D1*x.*sin(i)+3*x.*y.*cos(i)*sin(a)))./(4*pi*(x.^2+y.^2+D1.^2).^(5/2));
for k=1:nx
fprintf(fp,'%g ',Hax0(k));
end
fprintf(fp,'\n');
end
end
fclose(fp);
最后是将所有数据都输出至za_all.out以及hax_all.out文件中，其中有91（is从0-90°）*46（Deep从10-100m，间隔2m）行数据，每一行中都有nx（81，测线上的点数）个数据。
C#处理数据文件
首先将之前的za_all.out以及hax_all.out文件复制到WPF项目中，并且对其设置属性为“内容”，并“始终复制”。
将数据读取后分别存储在一个三维数组中，即倾角*深度*测点这样的结构。
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using System.IO;
using System.Resources;
using System.Windows;
using System.Windows.Resources;
namespace GeoMag_Viewer
{
class Values
{
private int iNum = 91;
private int jNum = 81;
public double[, ,] Za3D
{
get;
private set;
}
public double[, ,] Hax3D
{
get;
private set;
}
public Values()
{
StreamResourceInfo za3DInfo = Application.GetContentStream(new Uri("Resources/za_all.out", UriKind.Relative));
var za3DStreamReader = new StreamReader(za3DInfo.Stream);
string content = za3DStreamReader.ReadToEnd();
string[] lines = content.Split(new string[] { " \n" }, StringSplitOptions.None);
Za3D = new Double[91, 46, jNum];
for (int l = 0; l < 91; l++)
{
for (int i = 0; i < 46; i++)
{
int index = l * 46 + i;
string[] words = lines[index].Split(new string[] { " " }, StringSplitOptions.None);
for (int j = 0; j < jNum; j++)
{
Za3D[l, i, j] = Double.Parse(words[j]);
}
}
}
StreamResourceInfo hax3DInfo = Application.GetContentStream(new Uri("Resources/hax_all.out", UriKind.Relative));
var hax3DStreamReader = new StreamReader(hax3DInfo.Stream);
content = hax3DStreamReader.ReadToEnd();
lines = content.Split(new string[] { " \n" }, StringSplitOptions.None);
Hax3D = new Double[91, 46, jNum];
for (int l = 0; l < 91; l++)
{
for (int i = 0; i < 46; i++)
{
int index = l * 46 + i;
string[] words = lines[index].Split(new string[] { " " }, StringSplitOptions.None);
for (int j = 0; j < jNum; j++)
{
Hax3D[l, i, j] = Double.Parse(words[j]);
}
}
}
}
}
}使用OxyPlot组件
为了在WPF中绘制曲线图，通过搜索找到了OxyPlot这个开源库，可以方便地绘制我想要的图形。
创建WPF项目后，打开NUGet程序包管理器，搜索OxyPlot，安装OxyPlot for wpf。
Xaml代码如下：
<Window x:Class="GeoMag_Viewer.MainWindow"
xmlns="http://schemas.microsoft.com/winfx/2006/xaml/presentation"
xmlns:x="http://schemas.microsoft.com/winfx/2006/xaml"
xmlns:oxy="http://oxyplot.org/wpf"
xmlns:local="clr-namespace:GeoMag_Viewer"
Title="GeoMag Viewer" Height="350" Width="525">
<Grid Background="Gray">
<Grid.RowDefinitions>
<RowDefinition Height="*"/>
<RowDefinition Height="Auto"/>
<RowDefinition Height="Auto"/>
</Grid.RowDefinitions>
<Grid>
<Grid.ColumnDefinitions>
<ColumnDefinition Width="*"/>
<ColumnDefinition Width="Auto"/>
<ColumnDefinition Width="Auto"/>
</Grid.ColumnDefinitions>
<oxy:PlotView Grid.Column="0" Name="MyPlotView">
<oxy:PlotView.Series>
<oxy:LineSeries Name="ZaLine"/>
<oxy:LineSeries Name="HaxLine"/>
</oxy:PlotView.Series>
</oxy:PlotView>
<Image Grid.Column="0"
Name="SphereImage"
Source="Resources/ci.jpg"
HorizontalAlignment="Left"
Margin="150,0,0,0"
Width="150"
Height="150"
/>
<Slider Grid.Column="1" Name="DeepSlider"
Orientation="Vertical"
Margin="5,10,5,10"
ValueChanged="DeepSlider_ValueChanged"
Maximum="45"
Minimum="0"
TickFrequency="1"
TickPlacement="TopLeft"
/>
<TextBlock Grid.Column="2"
Name="DeepText"
FontSize="20"
FontWeight="Bold"
TextWrapping="Wrap"
VerticalAlignment="Center"
/>
</Grid>
<Slider Grid.Row="1"
Name="angleSlider"
ValueChanged="angleSlider_ValueChanged"
Margin="10,5,10,5"
Maximum="90"
Minimum="0"
TickFrequency="1"
TickPlacement="TopLeft"/>
<TextBlock Grid.Row="2"
Name="tickText"
FontSize="24"
FontWeight="Bold"
Margin="150,0,150,0"
HorizontalAlignment="Center"/>
</Grid>
</Window>
MainWindow.xaml.cs：
OxyPlot绘制曲线是是将Line.ItemSource设置为List，在这儿使用一个列表的二维数组分别保存Za、Hax异常。数组的第一个下标表示角度，第二个下标表示深度。
当用于控制角度、深度的Slider值变化时，分别调用angleSlider_ValueChanged和DeepSlider_ValueChanged来改变所要呈现的曲线，即
ZaLine.ItemsSource = Za3DLists[angleIndex, deepIndex];
HaxLine.ItemsSource = Hax3DLists[angleIndex, deepIndex];
这样，拖动Slider滑块，异常曲线就能跟着动态变化了，达到了演示效果。
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using System.Windows;
using System.Windows.Controls;
using System.Windows.Data;
using System.Windows.Documents;
using System.Windows.Input;
using System.Windows.Media;
using System.Windows.Media.Imaging;
using System.Windows.Navigation;
using System.Windows.Shapes;
using OxyPlot;
using OxyPlot.Series;
using OxyPlot.Wpf;
using Microsoft.Win32;
namespace GeoMag_Viewer
{
/// <summary>
/// MainWindow.xaml 的交互逻辑
/// </summary>
public partial class MainWindow : Window
{
public MainWindow()
{
InitializeComponent();
this.WindowState = WindowState.Maximized;
Values v = new Values();
MyPlotView.Title = "Za/Hax与磁化倾角的关系Demo";
var valueAxis = new LinearAxis() { Title = "磁异常/nT" };
MyPlotView.Axes.Add(valueAxis);
MyPlotView.LegendFontSize = 20;
ZaLine.Title = "Za异常";
HaxLine.Title = "Hax异常";
tickText.Text = "磁化倾角is=0";
DeepText.Text = "球体埋深";
double[, ,] Za3D = v.Za3D;
double[, ,] Hax3D = v.Hax3D;
int lNum = Za3D.GetLength(0);//角度数
int iNum = Za3D.GetLength(1);//深度数
int jNum = Za3D.GetLength(2);//单测线点数
Za3DLists = new List<DataPoint>[lNum, iNum];
Hax3DLists = new List<DataPoint>[lNum, iNum];
for (int l = 0; l < lNum; l++)
{
for (int i = 0; i < iNum; i++)
{
Za3DLists[l, i] = new List<DataPoint>();
Hax3DLists[l, i] = new List<DataPoint>();
for (int j = 0; j < jNum; j++)
{
int x = -200 + j * 5;
Za3DLists[l, i].Add(new DataPoint(x, Za3D[l, i, j]));
Hax3DLists[l, i].Add(new DataPoint(x, Hax3D[l, i, j]));
}
}
}
angleIndex = 0;
deepIndex = 0;
ZaLine.ItemsSource = Za3DLists[angleIndex, deepIndex];
HaxLine.ItemsSource = Hax3DLists[angleIndex, deepIndex];
width = 100;
height=100;
RotateTransform rt = new RotateTransform(0,width/2,height/2);
SphereImage.RenderTransform = rt;
}
double width;
double height;
private void angleSlider_ValueChanged(object sender, RoutedPropertyChangedEventArgs<double> e)
{
if (tickText == null)
{
return;
}
int angle = (int)e.NewValue;
int iAngle = angle;
angleIndex = iAngle;
tickText.Text = "磁化倾角is=" + angle;
ZaLine.ItemsSource = Za3DLists[angleIndex, deepIndex];
HaxLine.ItemsSource = Hax3DLists[angleIndex, deepIndex];
//旋转球体图像
//width / 2, height / 2 使绕图像中心旋转
RotateTransform rt = new RotateTransform(angle, width / 2, height / 2);
SphereImage.RenderTransform = rt;
}
private void DeepSlider_ValueChanged(object sender, RoutedPropertyChangedEventArgs<double> e)
{
if (DeepText == null)
{
return;
}
int iDeep = (int)e.NewValue;
int deep = 10 + 2 * iDeep;
deepIndex = iDeep;
DeepText.Text = "球体埋深\n" + deep+"m";
ZaLine.ItemsSource = Za3DLists[angleIndex, deepIndex];
HaxLine.ItemsSource = Hax3DLists[angleIndex, deepIndex];
}
//角度*深度
public List<DataPoint>[,] Za3DLists
{
get;
private set;
}
//角度*深度
public List<DataPoint>[,] Hax3DLists
{
get;
private set;
}
public int angleIndex
{
get;
private set;
}
public int deepIndex
{
get;
private set;
}
}
}另外，我还在程序中用了一个球体磁场图像，让它随着角度Slider的改变而旋转相应的角度，其实现也很简单，再次不再赘述。
版权声明：本文为博主原创文章，未经博主允许不得转载。