样条曲线。。。。。。。。。。

freshair · 发表于 2003-10-19 10:54:00

Private Sub Form_Load()
Const NC = 30, EE = 1E-30
Dim Y0(NC) '数组存放型值点的Y值
Dim Z0(NC) '数组存放型值点的Z值
Dim QY(NC), QZ(NC) '控制多边形上的点的坐标
Dim YTEM(NC), ZTEM(NC) '临时数组
Dim pathf As String '路径
Dim NY As Integer '型值点的个数
Dim YG As Single '给定的Y值
Dim ZQ As Single '待求的Z值
Dim ZG As Single '给定的Z值
Dim YQ As Single '待求的Y值
Dim I, J As Integer
pathf = App.Path
Open pathf & "\XING.TXT " For Input As #1 '已知型值点的输入文件
Open pathf & "\OUTPUT.TXT " For Append As #2
Input #1, NY
For I = 2 To NY + 1
Input #1, Y0(I), Z0(I)
'Print #2, Y0(I), Z0(I)
Next
Y0(1) = 0#
Z0(1) = 0#
Y0(NY + 2) = 0#
Z0(NY + 2) = 0#
'--------------------------------------------
Call FSQ(NY, Y0(), Z0(), QY(), QZ())
'-------------------------------------------
For I = 1 To 7
Input #1, YG
Call YTOZ(NY, YG, ZQ, Y0(), Z0(), QY(), QZ(), EE)
Print #2, YG, ",", ZQ
Next
'-------------------------------------------
For J = 1 To 10
Input #1, ZG
Call ZTOY(NY, ZG, YQ, Y0(), Z0(), QY(), QZ(), EE)
Print #2, YQ, ",", ZG
Next
Close #1
Close #2
End Sub
'反算由曲线上的点求控制多边形上的点
Sub FSQ(NY, Y0(), Z0(), QY(), QZ())
Const NC = 30
Dim AA(NC), BB(NC), CC(NC) '矩阵的系数
For I = 1 To NY + 2
Print #2, Y0(I), Z0(I)
Next
For I = 2 To NY + 1
AA(I) = 1#
BB(I) = 4#
CC(I) = 1#
Next
AA(NY + 2) = 6#
BB(1) = 6#
BB(NY + 2) = -6#
CC(1) = -6#
'For I = 1 To NY + 2
'Print #2, AA(I), BB(I), CC(I) '
'Next
For I = 1 To NY + 2
Y0(I) = Y0(I) * 6#
Z0(I) = Z0(I) * 6#
Next
Call TRIDAG(AA(), BB(), CC(), Y0(), QY(), (NY + 2)) '调用追赶法子程序
Call TRIDAG(AA(), BB(), CC(), Z0(), QZ(), (NY + 2)) '求控制多边形上的点坐标
'For I = 1 To (NY + 2)
'Print #2, QY(I), QZ(I)
'Next
For I = 1 To NY + 2
Y0(I) = Y0(I) / 6#
Z0(I) = Z0(I) / 6#
Next
End Sub
Sub YZT(YZG, YZ0(), DT, EE) '迭代求T
Dim XX1, XX2, XX, FF1, FF2, FF As Single
XX1 = 0#
XX2 = 1#
FF1 = ((-XX1 ^ 3 + 3 * XX1 ^ 2 - 3 * XX1 + 1) * YZ0(1) + (3 * XX1 ^ 3 - 6 * XX1 ^ 2 + 4) * YZ0(2) + (-3 * XX1 ^ 3 + 3 * XX1 ^ 2 + 3 * XX1 + 1) * YZ0(3) + XX1 ^ 3 * YZ0(4)) / 6 - YZG
FF2 = ((-XX2 ^ 3 + 3 * XX2 ^ 2 - 3 * XX2 + 1) * YZ0(1) + (3 * XX2 ^ 3 - 6 * XX2 ^ 2 + 4) * YZ0(2) + (-3 * XX2 ^ 3 + 3 * XX2 ^ 2 + 3 * XX2 + 1) * YZ0(3) + XX2 ^ 3 * YZ0(4)) / 6 - YZG
If (Sgn(FF2) <> Sgn(FF1)) Then
Do
XX = (XX1 + XX2) / 2#
FF = ((-XX ^ 3 + 3 * XX ^ 2 - 3 * XX + 1) * YZ0(1) + (3 * XX ^ 3 - 6 * XX ^ 2 + 4) * YZ0(2) + (-3 * XX ^ 3 + 3 * XX ^ 2 + 3 * XX + 1) * YZ0(3) + XX ^ 3 * YZ0(4)) / 6 - YZG
If FF = 0# Then
DT = XX
End If
If (Sgn(FF1) = Sgn(FF)) Then
XX1 = XX
FF1 = FF
End If
If (Sgn(FF2) = Sgn(FF)) Then
XX2 = XX
FF2 = FF
End If
Loop While (Abs(XX2 - XX1) > EE And (FF > EE))
End If
DT = (XX1 + XX2) / 2#
End Sub
Sub TRIDAG(A(), B(), C(), R(), U(), N) '解矩阵
NMAX = 100
Dim GAM(100)
If B(1) = 0# Then Exit Sub
BET = B(1)
U(1) = R(1) / BET
For J = 2 To N
GAM(J) = C(J - 1) / BET
BET = B(J) - A(J) * GAM(J)
If BET = 0# Then Exit Sub
U(J) = (R(J) - A(J) * U(J - 1)) / BET
Next J
For J = N - 1 To 1 Step -1
U(J) = U(J) - GAM(J + 1) * U(J + 1)
Next J
End Sub
'正算用矩阵法，给定Y求Z值
Sub YTOZ(NY, YG, ZQ, Y0(), Z0(), QY(), QZ(), EE)
Const NC = 30
Dim YZ1(4) '用于计算T值
Dim NI, II, J As Integer '迭待点所在的段
Dim XT As Single
Call LOCATE(Y0(), NY + 1, YG, J%)
NI = J
'Print #2, YG, NI '待求点在NI段上
'求T值
YZ1(1) = QY(NI - 1)
YZ1(2) = QY(NI)
YZ1(3) = QY(NI + 1)
YZ1(4) = QY(NI + 2)
Call YZT(YG, YZ1(), XT, EE)
'Print #2, T
ZQ = ((-XT ^ 3 + 3 * XT ^ 2 - 3 * XT + 1) * QZ(NI - 1) + (3 * XT ^ 3 - 6 * XT ^ 2 + 4) * QZ(NI) + (-3 * XT ^ 3 + 3 * XT ^ 2 + 3 * XT + 1) * QZ(NI + 1) + XT ^ 3 * QZ(NI + 2)) / 6
For II = 2 To NY + 1
If YG = Y0(II) Then
ZQ = Z0(II)
End If
Next
End Sub
Sub ZTOY(NY, ZG, YQ, Y0(), Z0(), QY(), QZ(), EE)
'正算用矩阵法，给定Z求Y值
Const NC = 30
Dim YZ2(4) '用于计算T值
Dim NI As Integer '迭待点所在的段数
Dim YT As Single
Call LOCATE(Z0(), NY + 1, ZG, J%)
NI = J
'Print #2, ZG, NI '待求点在NI段上
YZ2(1) = QZ(NI - 1)
YZ2(2) = QZ(NI)
YZ2(3) = QZ(NI + 1)
YZ2(4) = QZ(NI + 2)
Call YZT(ZG, YZ2(), YT, EE)
YQ = ((-YT ^ 3 + 3 * YT ^ 2 - 3 * YT + 1) * QY(NI - 1) + (3 * YT ^ 3 - 6 * YT ^ 2 + 4) * QY(NI) + (-3 * YT ^ 3 + 3 * YT ^ 2 + 3 * YT + 1) * QY(NI + 1) + YT ^ 3 * QY(NI + 2)) / 6
For II = 2 To NY + 1
If ZG = Z0(II) Then
YQ = Y0(II)
End If
Next
End Sub
Sub LOCATE(XX(), N, X, J%) '判断待求点所在段号
JL = 0
JU = N + 1
10 If JU - JL > 1 Then
JM = (JU + JL) / 2
If XX(N) > XX(1) Eqv X > XX(JM) Then
JL = JM
Else
JU = JM
End If
GoTo 10
End If
J% = JL
End Sub
===================================
比如输入文件xing.text如下
13
2.9193 0.12086
4.500 0.3233
5.117 0.50
6.000 0.8965
6.8289 1.50
7.6308 2.50
8.1327 3.50
8.4930 4.50
8.6820 5.20
8.7508 5.50
8.9382 6.50
9.0758 7.50
9.1773 8.50
3
4
5
6
7
8
9
0.2
0.8
1
2
3
4
5
6
7
8

Private Sub Form_Load() 
Const NC = 30, EE = 1E-30 
Dim Y0(NC) '数组存放型值点的Y值 
Dim Z0(NC) '数组存放型值点的Z值 
Dim QY(NC), QZ(NC) '控制多边形上的点的坐标 
Dim YTEM(NC), ZTEM(NC) '临时数组 
Dim pathf As String '路径 
Dim NY As Integer '型值点的个数 
Dim YG As Single '给定的Y值 
Dim ZQ As Single '待求的Z值 
Dim ZG As Single '给定的Z值 
Dim YQ As Single '待求的Y值 
Dim I, J As Integer 
 
 
pathf = App.Path 
Open pathf & "\XING.TXT " For Input As #1 '已知型值点的输入文件 
Open pathf & "\OUTPUT.TXT " For Append As #2 
Input #1, NY 
For I = 2 To NY + 1 
Input #1, Y0(I), Z0(I) 
'Print #2, Y0(I), Z0(I) 
Next 
Y0(1) = 0# 
Z0(1) = 0# 
Y0(NY + 2) = 0# 
Z0(NY + 2) = 0# 
'-------------------------------------------- 
Call FSQ(NY, Y0(), Z0(), QY(), QZ()) 
'------------------------------------------- 
For I = 1 To 7 
Input #1, YG 
Call YTOZ(NY, YG, ZQ, Y0(), Z0(), QY(), QZ(), EE) 
Print #2, YG, ",", ZQ 
Next 
'------------------------------------------- 
For J = 1 To 10 
Input #1, ZG 
Call ZTOY(NY, ZG, YQ, Y0(), Z0(), QY(), QZ(), EE) 
Print #2, YQ, ",", ZG 
Next 
Close #1 
Close #2 
End Sub 
 
'反算由曲线上的点求控制多边形上的点 
Sub FSQ(NY, Y0(), Z0(), QY(), QZ()) 
Const NC = 30 
Dim AA(NC), BB(NC), CC(NC) '矩阵的系数 
 
For I = 1 To NY + 2 
Print #2, Y0(I), Z0(I) 
Next 
For I = 2 To NY + 1 
AA(I) = 1# 
BB(I) = 4# 
CC(I) = 1# 
Next 
AA(NY + 2) = 6# 
BB(1) = 6# 
BB(NY + 2) = -6# 
CC(1) = -6# 
'For I = 1 To NY + 2 
'Print #2, AA(I), BB(I), CC(I) ' 
'Next 
For I = 1 To NY + 2 
Y0(I) = Y0(I) * 6# 
Z0(I) = Z0(I) * 6# 
Next 
 
Call TRIDAG(AA(), BB(), CC(), Y0(), QY(), (NY + 2))  '调用追赶法子程序 
Call TRIDAG(AA(), BB(), CC(), Z0(), QZ(), (NY + 2))  '求控制多边形上的点坐标 
'For I = 1 To (NY + 2) 
'Print #2, QY(I), QZ(I) 
'Next 
For I = 1 To NY + 2 
Y0(I) = Y0(I) / 6# 
Z0(I) = Z0(I) / 6# 
Next 
End Sub 
 
Sub YZT(YZG, YZ0(), DT, EE) '迭代求T 
Dim XX1, XX2, XX, FF1, FF2, FF As Single 
XX1 = 0# 
XX2 = 1# 
 
FF1 = ((-XX1 ^ 3 + 3 * XX1 ^ 2 - 3 * XX1 + 1) * YZ0(1) + (3 * XX1 ^ 3 - 6 * XX1 ^ 2 + 4) * YZ0(2) + (-3 * XX1 ^ 3 + 3 * XX1 ^ 2 + 3 * XX1 + 1) * YZ0(3) + XX1 ^ 3 * YZ0(4)) / 6 - YZG 
FF2 = ((-XX2 ^ 3 + 3 * XX2 ^ 2 - 3 * XX2 + 1) * YZ0(1) + (3 * XX2 ^ 3 - 6 * XX2 ^ 2 + 4) * YZ0(2) + (-3 * XX2 ^ 3 + 3 * XX2 ^ 2 + 3 * XX2 + 1) * YZ0(3) + XX2 ^ 3 * YZ0(4)) / 6 - YZG 
 
If (Sgn(FF2) <> Sgn(FF1)) Then 
Do 
 
XX = (XX1 + XX2) / 2# 
FF = ((-XX ^ 3 + 3 * XX ^ 2 - 3 * XX + 1) * YZ0(1) + (3 * XX ^ 3 - 6 * XX ^ 2 + 4) * YZ0(2) + (-3 * XX ^ 3 + 3 * XX ^ 2 + 3 * XX + 1) * YZ0(3) + XX ^ 3 * YZ0(4)) / 6 - YZG 
If FF = 0# Then 
 DT = XX 
 End If 
If (Sgn(FF1) = Sgn(FF)) Then 
 XX1 = XX 
 FF1 = FF 
 End If 
If (Sgn(FF2) = Sgn(FF)) Then 
XX2 = XX 
FF2 = FF 
End If 
Loop While (Abs(XX2 - XX1) > EE And (FF > EE)) 
 
End If 
 
DT = (XX1 + XX2) / 2# 
End Sub 
Sub TRIDAG(A(), B(), C(), R(), U(), N) '解矩阵 
    NMAX = 100 
    Dim GAM(100) 
    If B(1) = 0# Then Exit Sub 
    BET = B(1) 
    U(1) = R(1) / BET 
    For J = 2 To N 
        GAM(J) = C(J - 1) / BET 
        BET = B(J) - A(J) * GAM(J) 
        If BET = 0# Then Exit Sub 
        U(J) = (R(J) - A(J) * U(J - 1)) / BET 
    Next J 
    For J = N - 1 To 1 Step -1 
        U(J) = U(J) - GAM(J + 1) * U(J + 1) 
    Next J 
End Sub 
'正算用矩阵法，给定Y求Z值 
Sub YTOZ(NY, YG, ZQ, Y0(), Z0(), QY(), QZ(), EE) 
Const NC = 30 
Dim YZ1(4) '用于计算T值 
Dim NI, II, J As Integer '迭待点所在的段 
Dim XT As Single 
 
 
 
Call LOCATE(Y0(), NY + 1, YG, J%) 
 
NI = J 
 
'Print #2, YG, NI '待求点在NI段上 
 
 
 
'求T值 
YZ1(1) = QY(NI - 1) 
YZ1(2) = QY(NI) 
YZ1(3) = QY(NI + 1) 
YZ1(4) = QY(NI + 2) 
Call YZT(YG, YZ1(), XT, EE) 
 
'Print #2, T 
ZQ = ((-XT ^ 3 + 3 * XT ^ 2 - 3 * XT + 1) * QZ(NI - 1) + (3 * XT ^ 3 - 6 * XT ^ 2 + 4) * QZ(NI) + (-3 * XT ^ 3 + 3 * XT ^ 2 + 3 * XT + 1) * QZ(NI + 1) + XT ^ 3 * QZ(NI + 2)) / 6 
For II = 2 To NY + 1 
 If YG = Y0(II) Then 
ZQ = Z0(II) 
 End If 
Next 
End Sub 
Sub ZTOY(NY, ZG, YQ, Y0(), Z0(), QY(), QZ(), EE) 
'正算用矩阵法，给定Z求Y值 
Const NC = 30 
Dim YZ2(4) '用于计算T值 
Dim NI As Integer '迭待点所在的段数 
Dim YT As Single 
 
 
Call LOCATE(Z0(), NY + 1, ZG, J%) 
NI = J 
'Print #2, ZG, NI '待求点在NI段上 
 
YZ2(1) = QZ(NI - 1) 
YZ2(2) = QZ(NI) 
YZ2(3) = QZ(NI + 1) 
YZ2(4) = QZ(NI + 2) 
Call YZT(ZG, YZ2(), YT, EE) 
 
YQ = ((-YT ^ 3 + 3 * YT ^ 2 - 3 * YT + 1) * QY(NI - 1) + (3 * YT ^ 3 - 6 * YT ^ 2 + 4) * QY(NI) + (-3 * YT ^ 3 + 3 * YT ^ 2 + 3 * YT + 1) * QY(NI + 1) + YT ^ 3 * QY(NI + 2)) / 6 
For II = 2 To NY + 1 
 If ZG = Z0(II) Then 
YQ = Y0(II) 
 End If 
Next 
End Sub 
Sub LOCATE(XX(), N, X, J%) '判断待求点所在段号 
    JL = 0 
    JU = N + 1 
10  If JU - JL > 1 Then 
        JM = (JU + JL) / 2 
        If XX(N) > XX(1) Eqv X > XX(JM) Then 
            JL = JM 
        Else 
            JU = JM 
        End If 
        GoTo 10 
    End If 
    J% = JL 
End Sub 
 
=================================== 
比如输入文件xing.text如下 
 
13 
2.9193 0.12086 
4.500 0.3233 
5.117 0.50 
6.000 0.8965 
6.8289 1.50 
7.6308 2.50 
8.1327 3.50 
8.4930 4.50 
8.6820 5.20 
8.7508 5.50 
8.9382 6.50 
9.0758 7.50 
9.1773 8.50 
3 
4 
5 
6 
7 
8 
9 
0.2 
0.8 
1 
2 
3 
4 
5 
6 
7 
8

gzy · 发表于 2003-10-19 18:20:00

又有人往上贴了，这好象是继明总之后的关于VBA的第一个。
不过freshairMM，这叫源码共享，你这好象还不够“源”啊，现在这里人气不旺，你就全部拿出来吧。

freshair · 发表于 2003-10-19 21:55:00

呵呵，已经比较源了，我没有保留啊。今晚来不及了，如有问题，我明天再补上。

asahix · 发表于 2003-11-21 13:31:00

[求助]freshair好象和我是同行，我是做船厂船舶设计的，请问有没有曲线拟合的程序？
就是给出一组数据把需要的比如三次多项式或者乘幂的曲线得出来。

无痕 · 发表于 2003-12-22 06:50:00

楼主的程序用在什么地方呢?

BDYCAD · 发表于 2005-11-1 14:00:00

要是可以有象CorelDRAW 畫的曲線多好啊!

账号		自动登录	找回密码
密码			注册