C#地震反应谱
好久没有发东西了，把一年多前写的C#的地震反应谱代码贴一下，供大家参考。
using System;
using System.Collections.Generic;
using System.ComponentModel;
using System.Data;
using System.Drawing;
using System.Linq;
using System.Text;
using System.Windows.Forms;
using ZedGraph;
using System.IO;
namespace 地震反应谱计算
{
public partial class Form1 : Form
{
PointPairList list0 = new PointPairList();
PointPairList list1 = new PointPairList();
PointPairList list2 = new PointPairList();
PointPairList list3 = new PointPairList();
int signumb;
double[] tim = new double[10000];
double[] vibsign = new double[10000];
double[] tempdis = new double[10000];
double[] tempvel = new double[10000];
double[] tempacc = new double[10000];
double[] accespe = new double[2000];
double[] velospe = new double[2000];
double[] dispspe = new double[2000];
double[] freoper = new double[2000];
double maxperi, minperi, delperi;
double[] z = new double[3];
double pi = 3.14159;
int sampfre, i, j, spenumb;
string filename = "";
double max_y = 0, min_y = 0;
public Form1()
{
InitializeComponent();
}
private void Form1_Load(object sender, EventArgs e)
{
using (StreamReader sr = new StreamReader(new FileStream("initial.dat", FileMode.Open, FileAccess.Read)))
{
i = 0;
string tmp = "";
string[] ini = new string[9];
while (sr.Peek() > 0)
{
string[] t = new string[2];
tmp = sr.ReadLine();
t = tmp.Split('\t');
ini[i] =t[0];
i++;
}
filename = ini[0];
sampfre = Convert.ToInt16(ini[1]);
signumb = Convert.ToInt16(ini[2]);
minperi = Convert.ToDouble(ini[3]);
maxperi = Convert.ToDouble(ini[4]);
delperi = Convert.ToDouble(ini[5]);
z[0] = Convert.ToDouble(ini[6]);
z[1] = Convert.ToDouble(ini[7]);
z[2] = Convert.ToDouble(ini[8]);
}
spenumb = Convert.ToInt16((maxperi - minperi) / delperi + 1);
for (i = 0; i < spenumb; i++)
{
freoper[i] = minperi + delperi*i;
}
using (StreamReader sr = new StreamReader(new FileStream(filename, FileMode.Open, FileAccess.Read)))
{
i = 0;
string tmp = "";
while (sr.Peek() > 0)
{
string[] t = new string[2];
tmp = sr.ReadLine();
t = tmp.Split('\t');
vibsign[i] = Convert.ToDouble(t[1]);
tim[i] = Convert.ToDouble(t[0]);
list0.Add(tim[i], vibsign[i]);
max_y = Math.Max(vibsign[i], max_y);
min_y = Math.Min(vibsign[i], min_y);
i++;
}
}
GraphPane myPane = zg1.GraphPane;
myPane.XAxis.Scale.Min = 0;
myPane.XAxis.Scale.Max =signumb/sampfre;
myPane.YAxis.Scale.Min = min_y;
myPane.YAxis.Scale.Max = max_y;
myPane.Title.Text = "时程曲线";
myPane.XAxis.Title.Text = "时间/s";
myPane.YAxis.Title.Text = "加速度/(m/s\u00b2)";
LineItem myCurve0 = myPane.AddCurve("时程曲线", list0, Color.Red, SymbolType.None);
myCurve0.Symbol.Fill = new Fill(Color.White);
zg1.Refresh();
}
private void button1_Click(object sender, EventArgs e)
{
zg1.SaveAs();
}
private void radioButton1_CheckedChanged(object sender, EventArgs e)
{
GraphPane myPane = zg1.GraphPane;
myPane.XAxis.Scale.Min = 0;
myPane.XAxis.Scale.Max = signumb / sampfre;
myPane.YAxis.Scale.Min = min_y;
myPane.YAxis.Scale.Max = max_y;
myPane.CurveList.Clear();
list1.Clear();
list2.Clear();
list3.Clear();
myPane.Title.Text = "时程曲线";
myPane.XAxis.Title.Text = "时间/s";
myPane.YAxis.Title.Text = "加速度/(m/s\u00b2)";
LineItem myCurve0 = myPane.AddCurve("时程曲线", list0, Color.Red, SymbolType.None);
myCurve0.Symbol.Fill = new Fill(Color.White);
zg1.Refresh();
}
private void button2_Click(object sender, EventArgs e)
{
zg1.DoPrintPreview();
}
private void button3_Click(object sender, EventArgs e)
{
zg1.DoPageSetup();
}
private void button4_Click(object sender, EventArgs e)
{
zg1.DoPrint();
}
private void radioButton2_CheckedChanged(object sender, EventArgs e)
{
double maxy = 0;
GraphPane myPane = zg1.GraphPane;
myPane.CurveList.Clear();
list1.Clear();
list2.Clear();
list3.Clear();
myPane.Title.Text = "加速度反应谱";
myPane.XAxis.Title.Text = "周期/s";
myPane.YAxis.Title.Text = "加速度/(m/s\u00b2)";
int type;
for (type = 0; type < 3; type++)
{
for (i = 0; i < spenumb; i++)
{
dispspe[i] = 0.0;
velospe[i] = 0.0;
accespe[i] = 0.0;
}
for (i = 0; i < spenumb; i++)
{
double damprat = z[type];
double angfreq = 2 * pi / freoper[i]; ;
double deltime = 1.0 / sampfre;
double a = angfreq * Math.Sqrt(1 - (damprat * damprat)) * deltime;
double b = damprat * angfreq * deltime;
double d = Math.Sqrt(1 - (damprat * damprat));
double ex = Math.Exp(-b);
double a11 = ex * (damprat * Math.Sin(a) / d + Math.Cos(a));
double a12 = ex * Math.Sin(a) / angfreq / d;
double a21 = ex * (-angfreq) * Math.Sin(a) / d;
double a22 = ex * (-damprat * Math.Sin(a) / d + Math.Cos(a));
double b11 = ex * ((2 * (damprat * damprat) - 1 + b) * Math.Sin(a) / (angfreq * angfreq) / a + (2 * damprat + angfreq * deltime) * Math.Cos(a) / (angfreq * angfreq * angfreq) / deltime) - 2 * damprat / (angfreq * angfreq * angfreq) / deltime;
double b12 = ex * ((1 - 2 * (damprat * damprat)) * Math.Sin(a) / (angfreq * angfreq) / a - 2 * damprat * Math.Cos(a) / (angfreq * angfreq * angfreq) / deltime) - 1 / (angfreq * angfreq) + 2 * damprat / (angfreq * angfreq * angfreq) / deltime;
double b21 = ex * ((-damprat - angfreq * deltime) * Math.Sin(a) / angfreq / a - Math.Cos(a) / (angfreq * angfreq) / deltime) + 1 / (angfreq * angfreq) / deltime;
double b22 = ex * (damprat * Math.Sin(a) / angfreq / a + Math.Cos(a) / (angfreq * angfreq) / deltime) - 1 / (angfreq * angfreq) / deltime;
tempdis[0] = 0;
tempvel[0] = 0;
tempacc[0] = 0;
for (j = 1; j < signumb; j++)
{
tempdis[j] = a11 * tempdis[j - 1] + a12 * tempvel[j - 1] + b11 * vibsign[j - 1] + b12 * vibsign[j];
tempvel[j] = a21 * tempdis[j - 1] + a22 * tempvel[j - 1] + b21 * vibsign[j - 1] + b22 * vibsign[j];
tempacc[j] = -(2 * damprat * angfreq * tempvel[j] + angfreq * angfreq * tempdis[j]);
dispspe[i] = Math.Max(dispspe[i], Math.Abs(tempdis[j]));
velospe[i] = Math.Max(velospe[i], Math.Abs(tempvel[j]));
accespe[i] = Math.Max(accespe[i], Math.Abs(tempacc[j]));
maxy = Math.Max(accespe[i],maxy);
}
if (type == 0)
{
list1.Add(delperi * i, accespe[i]);
}
else if (type == 1)
{
list2.Add(delperi * i, accespe[i]);
}
else if (type == 2)
{
list3.Add(delperi * i, accespe[i]);
}
}
}
myPane.XAxis.Scale.Min = 0;
myPane.XAxis.Scale.Max = maxperi;
myPane.YAxis.Scale.Min = 0;
myPane.YAxis.Scale.Max = maxy;
LineItem myCurve1 = myPane.AddCurve("ξ=" + z[0].ToString(), list1, Color.Red, SymbolType.None);
myCurve1.Symbol.Fill = new Fill(Color.White);
LineItem myCurve2 = myPane.AddCurve("ξ=" + z[1].ToString(), list2, Color.Green, SymbolType.None);
myCurve2.Symbol.Fill = new Fill(Color.White);
LineItem myCurve3 = myPane.AddCurve("ξ=" + z[2].ToString(), list3, Color.Blue, SymbolType.None);
myCurve3.Symbol.Fill = new Fill(Color.White);
zg1.Refresh();
}
private void radioButton3_CheckedChanged(object sender, EventArgs e)
{
double maxy = 0;
GraphPane myPane = zg1.GraphPane;
myPane.CurveList.Clear();
list1.Clear();
list2.Clear();
list3.Clear();
myPane.Title.Text = "速度反应谱";
myPane.XAxis.Title.Text = "周期/s";
myPane.YAxis.Title.Text = "速度/(cm/s)";
int type;
for (type = 0; type < 3; type++)
{
for (i = 0; i < spenumb; i++)
{
dispspe[i] = 0.0;
velospe[i] = 0.0;
accespe[i] = 0.0;
}
for (i = 0; i < spenumb; i++)
{
double damprat = z[type];
double angfreq = 2 * pi / freoper[i]; ;
double deltime = 1.0 / sampfre;
double a = angfreq * Math.Sqrt(1 - (damprat * damprat)) * deltime;
double b = damprat * angfreq * deltime;
double d = Math.Sqrt(1 - (damprat * damprat));
double ex = Math.Exp(-b);
double a11 = ex * (damprat * Math.Sin(a) / d + Math.Cos(a));
double a12 = ex * Math.Sin(a) / angfreq / d;
double a21 = ex * (-angfreq) * Math.Sin(a) / d;
double a22 = ex * (-damprat * Math.Sin(a) / d + Math.Cos(a));
double b11 = ex * ((2 * (damprat * damprat) - 1 + b) * Math.Sin(a) / (angfreq * angfreq) / a + (2 * damprat + angfreq * deltime) * Math.Cos(a) / (angfreq * angfreq * angfreq) / deltime) - 2 * damprat / (angfreq * angfreq * angfreq) / deltime;
double b12 = ex * ((1 - 2 * (damprat * damprat)) * Math.Sin(a) / (angfreq * angfreq) / a - 2 * damprat * Math.Cos(a) / (angfreq * angfreq * angfreq) / deltime) - 1 / (angfreq * angfreq) + 2 * damprat / (angfreq * angfreq * angfreq) / deltime;
double b21 = ex * ((-damprat - angfreq * deltime) * Math.Sin(a) / angfreq / a - Math.Cos(a) / (angfreq * angfreq) / deltime) + 1 / (angfreq * angfreq) / deltime;
double b22 = ex * (damprat * Math.Sin(a) / angfreq / a + Math.Cos(a) / (angfreq * angfreq) / deltime) - 1 / (angfreq * angfreq) / deltime;
tempdis[0] = 0;
tempvel[0] = 0;
tempacc[0] = 0;
for (j = 1; j < signumb; j++)
{
tempdis[j] = a11 * tempdis[j - 1] + a12 * tempvel[j - 1] + b11 * vibsign[j - 1] + b12 * vibsign[j];
tempvel[j] = a21 * tempdis[j - 1] + a22 * tempvel[j - 1] + b21 * vibsign[j - 1] + b22 * vibsign[j];
tempacc[j] = -(2 * damprat * angfreq * tempvel[j] + angfreq * angfreq * tempdis[j]);
dispspe[i] = Math.Max(dispspe[i], Math.Abs(tempdis[j]));
velospe[i] = Math.Max(velospe[i], Math.Abs(tempvel[j]));
accespe[i] = Math.Max(accespe[i], Math.Abs(tempacc[j]));
maxy = Math.Max(velospe[i], maxy);
}
if (type == 0)
{
list1.Add(delperi * i, velospe[i]*100);
}
else if (type == 1)
{
list2.Add(delperi * i, velospe[i] * 100);
}
else if (type == 2)
{
list3.Add(delperi * i, velospe[i] * 100);
}
}
}
myPane.XAxis.Scale.Min = 0;
myPane.XAxis.Scale.Max = maxperi;
myPane.YAxis.Scale.Min = 0;
myPane.YAxis.Scale.Max = maxy * 100;
LineItem myCurve1 = myPane.AddCurve("ξ=" + z[0].ToString(), list1, Color.Red, SymbolType.None);
myCurve1.Symbol.Fill = new Fill(Color.White);
LineItem myCurve2 = myPane.AddCurve("ξ=" + z[1].ToString(), list2, Color.Green, SymbolType.None);
myCurve2.Symbol.Fill = new Fill(Color.White);
LineItem myCurve3 = myPane.AddCurve("ξ=" + z[2].ToString(), list3, Color.Blue, SymbolType.None);
myCurve3.Symbol.Fill = new Fill(Color.White);
zg1.Refresh();
}
private void radioButton4_CheckedChanged(object sender, EventArgs e)
{
double maxy = 0;
GraphPane myPane = zg1.GraphPane;
myPane.CurveList.Clear();
list1.Clear();
list2.Clear();
list3.Clear();
myPane.Title.Text = "位移反应谱";
myPane.XAxis.Title.Text = "周期/s";
myPane.YAxis.Title.Text = "位移/(cm)";
int type;
for (type = 0; type < 3; type++)
{
for (i = 0; i < spenumb; i++)
{
dispspe[i] = 0.0;
velospe[i] = 0.0;
accespe[i] = 0.0;
}
for (i = 0; i < spenumb; i++)
{
double damprat = z[type];
double angfreq = 2 * pi / freoper[i]; ;
double deltime = 1.0 / sampfre;
double a = angfreq * Math.Sqrt(1 - (damprat * damprat)) * deltime;
double b = damprat * angfreq * deltime;
double d = Math.Sqrt(1 - (damprat * damprat));
double ex = Math.Exp(-b);
double a11 = ex * (damprat * Math.Sin(a) / d + Math.Cos(a));
double a12 = ex * Math.Sin(a) / angfreq / d;
double a21 = ex * (-angfreq) * Math.Sin(a) / d;
double a22 = ex * (-damprat * Math.Sin(a) / d + Math.Cos(a));
double b11 = ex * ((2 * (damprat * damprat) - 1 + b) * Math.Sin(a) / (angfreq * angfreq) / a + (2 * damprat + angfreq * deltime) * Math.Cos(a) / (angfreq * angfreq * angfreq) / deltime) - 2 * damprat / (angfreq * angfreq * angfreq) / deltime;
double b12 = ex * ((1 - 2 * (damprat * damprat)) * Math.Sin(a) / (angfreq * angfreq) / a - 2 * damprat * Math.Cos(a) / (angfreq * angfreq * angfreq) / deltime) - 1 / (angfreq * angfreq) + 2 * damprat / (angfreq * angfreq * angfreq) / deltime;
double b21 = ex * ((-damprat - angfreq * deltime) * Math.Sin(a) / angfreq / a - Math.Cos(a) / (angfreq * angfreq) / deltime) + 1 / (angfreq * angfreq) / deltime;
double b22 = ex * (damprat * Math.Sin(a) / angfreq / a + Math.Cos(a) / (angfreq * angfreq) / deltime) - 1 / (angfreq * angfreq) / deltime;
tempdis[0] = 0;
tempvel[0] = 0;
tempacc[0] = 0;
for (j = 1; j < signumb; j++)
{
tempdis[j] = a11 * tempdis[j - 1] + a12 * tempvel[j - 1] + b11 * vibsign[j - 1] + b12 * vibsign[j];
tempvel[j] = a21 * tempdis[j - 1] + a22 * tempvel[j - 1] + b21 * vibsign[j - 1] + b22 * vibsign[j];
tempacc[j] = -(2 * damprat * angfreq * tempvel[j] + angfreq * angfreq * tempdis[j]);
dispspe[i] = Math.Max(dispspe[i], Math.Abs(tempdis[j]));
velospe[i] = Math.Max(velospe[i], Math.Abs(tempvel[j]));
accespe[i] = Math.Max(accespe[i], Math.Abs(tempacc[j]));
maxy = Math.Max(dispspe[i], maxy);
}
if (type == 0)
{
list1.Add(delperi * i, dispspe[i]*100);
}
else if (type == 1)
{
list2.Add(delperi * i, dispspe[i] * 100);
}
else if (type == 2)
{
list3.Add(delperi * i, dispspe[i] * 100);
}
}
}
myPane.XAxis.Scale.Min = 0;
myPane.XAxis.Scale.Max = maxperi;
myPane.YAxis.Scale.Min = 0;
myPane.YAxis.Scale.Max = maxy*100;
LineItem myCurve1 = myPane.AddCurve("ξ=" + z[0].ToString(), list1, Color.Red, SymbolType.None);
myCurve1.Symbol.Fill = new Fill(Color.White);
LineItem myCurve2 = myPane.AddCurve("ξ=" + z[1].ToString(), list2, Color.Green, SymbolType.None);
myCurve2.Symbol.Fill = new Fill(Color.White);
LineItem myCurve3 = myPane.AddCurve("ξ=" + z[2].ToString(), list3, Color.Blue, SymbolType.None);
myCurve3.Symbol.Fill = new Fill(Color.White);
zg1.Refresh();
}
}
}
程序用了一个图形插件ZedGraph.dll
设置文件为initial.dat
input.dat %输入文件名 100 %采样频率 1500 %采样点数 0.01 %反应谱计算周期最小值 6 %反应谱计算周期最大值 0.01 %反应谱计算周期间隔 0.01 %计算阻尼1 0.05 %计算阻尼2 0.1 %计算阻尼3
地震波文件为input.dat，此处省略
