Public Type Complex x As Double y As Double End Type
Function ToComplex(x As Double, y As Double) As Complex ToComplex.x = x ToComplex.y = y End Function
Function Cadd(z1 As Complex, z2 As Complex) As Complex Cadd.x = z1.x + z2.x Cadd.y = z1.y + z2.y End Function
Function Csub(z1 As Complex, z2 As Complex) As Complex Csub.x = z1.x - z2.x Csub.y = z1.y - z2.y End Function
Function Cmul(z1 As Complex, z2 As Complex) As Complex Cmul.x = z1.x * z2.x - z1.y * z2.y Cmul.y = z1.y * z2.x + z1.x * z2.y End Function
Function Cdiv(z1 As Complex, z2 As Complex) As Complex Dim c As Double c = z2.x * z2.x + z2.y * z2.y Cdiv.x = (z1.x * z2.x + z1.y * z2.y) / c Cdiv.y = (z1.y * z2.x - z1.x * z2.y) / c
End Function
Function Conj(z As Complex) As Complex Conj.x = z.x comj.y = -z.y End Function
Function Rez(z As Complex) As Double Rez = z.x End Function
Function Imz(z As Complex) As Double Imz = z.y End Function
Function Cabs(z As Complex) As Double Cabs = Sqr(z.x * z.x + z.y * z.y) End Function
pi = 3.14159265358979 If x = 0# And y = 0# Then atan2 = 0# Exit Function End If
If x = 0# Then If y > 0 Then atan2 = pi Exit Function Else atan2 = -pi Exit Function End If End If
If x > 0 Then atan2 = Atn(y / x) Else If y > 0 Then atan2 = pi + Atn(y / x) Else atan2 = -pi + Atn(y / x) End If End If
End Function
Function Cexp(z As Complex) As Complex Cexp.x = Exp(z.x) * Cos(z.y) Cexp.y = Exp(z.x) * Sin(z.y) End Function
Function Clog(z As Complex) As Complex Dim R As Double R = z.x * z.x + z.y * z.y If R <> 0 Then Clog.x = 0.5 * Log(R) Clog.y = atan2(z.x, z.y) Else Clog.x = 0# Clog.y = 0# End If End Function
Function Cpow(z1 As Complex, z2 As Complex) As Complex Cpow = Cexp(Cmul(Clog(z1), z2)) End Function
Function ccosh(z As Complex) As Complex ccosh.x = (Exp(z.x) * Cos(z.y) + Exp(-z.x) * Cos(-z.y)) / 2 ccosh.y = (Exp(z.x) * Sin(z.y) + Exp(-z.x) * Sin(-z.y)) / 2 End Function
Function csinh(z As Complex) As Complex csinh.x = (Exp(z.x) * Cos(z.y) - Exp(-z.x) * Cos(-z.y)) / 2 csinh.y = (Exp(z.x) * Sin(z.y) - Exp(-z.x) * Sin(-z.y)) / 2 End Function
Rem 数値的逆ラプラス変換filt07(t) Hosonoのアルゴリズム-------------- Rem atstp:分母 btstp:分子---ラプラス関数値の算出に必要なパラメータをセル参照 Rem 作例 F(s) = ht(8次多項式)/gt(8次多項式) : 有理関数 Public Function filt07(time As Double, atstp As Variant, btstp As Variant, nts As Integer, p As Double) As Double Dim ss(10000) As Complex Dim x1L(10000) As Complex Dim res(10000) As Double Dim aa As Double Dim n As Double Dim cnt As Double
Rem パラメータ関連・複素数変数 Dim resigma As Double Dim ared(20, 0) As Double Dim bred(20, 0) As Double Dim pi As Double Dim a01(20) As Complex Dim b01(20) As Complex Dim ht(200) As Complex Dim gt(200) As Complex Dim ht01 As Complex Dim ht02 As Complex Dim ht03 As Complex Dim ht04 As Complex Dim gt01 As Complex Dim gt02 As Complex Dim gt03 As Complex Dim gt04 As Complex
Dim ap(500) As Double Dim cp(500) As Double
Rem 各種定数 pi = 3.14159265358979 Rem オイラー変換係数 cp(0) = 1 ap(p + 1) = 0 For q = 1 To p cp(q) = cp(q - 1) * (p + 1 - q) / q Next q For q = p To 0 Step -1 ap(q) = ap(q + 1) + cp(q) Next q
Rem パラメータの取り込み・複素数変換 Rem 指定範囲から順番にデータ用セル値を取り出す(次数の高い順a0-a8,b0-b8) For Each myCell In atstp ared(cnt, 0) = myCell.Value a01(cnt) = ToComplex(ared(cnt, 0), 0) cnt = cnt + 1 Next cnt = 0 For Each myCell In btstp bred(cnt, 0) = myCell.Value b01(cnt) = ToComplex(bred(cnt, 0), 0) cnt = cnt + 1 Next
Rem nts項分加算 resigma = 0 For n = 1 To nts resigma = resigma + res(n) Next
Rem オイラー変換-p項分追加 For i = 1 To p resigma = resigma + res(nts + i) * ap(i) / 2 ^ (p) Next Rem 変換値出力 filt07 = resigma * Exp(aa) / time Else Rem t=0の関数値:初期値の出力 filt07 = 0 End If End Function
Function csqr(z As Complex) As Complex Dim xx As Double Dim yy As Double xx = z.x yy = z.y R = Sqr((xx) ^ 2 + (yy) ^ 2) csqr.x = Sqr(xx + R) / Sqr(2) csqr.y = Sqr(R - xx) / Sqr(2) If z.y < 0 Then csqr.y = -csqr.y End Function