実行する際は、線形状（ベジェ曲線)を選択し、形状に合った条件を引数に代入してください。

スクリプト内に記述　詳細は上記のルーチン名をクリックして下さい。

同上

Option Explicit
XShade.Message"====================start====================="
Call Main()
XShade.Message"====================end====================="

'================================================================================
'==  曲線(ベジェ曲線）に関する共通 main
'================================================================================

'----------------ベジェ曲線用座標ﾃﾞｰﾀの取得
sub Main()
  Dim NA,Vec3,PtNo
  Dim Rou,a,b,c,ta,tb,n,Ltarget
  Dim eVec,Pos,L,nd
  Dim tt,NoOfT,i,DatXYZ,Direction

  'ｱﾝｶｰﾎﾟｲﾝﾄの数を取得
  NA=XShade.NumberOfAnchorPoints()
  'ｱﾝｶｰポイントの座標値を取得
  Vec3=GetVec3(NA)

'--------------------ベジェ曲線でt=taの地点の座標を計算する。
'    (ｱﾝｶｰﾎﾟｲﾝﾄ番号PtNo〜PtNo+1の部分が対象）
'  ta=0.32: PtNo=0
'  Pos=calcPos(Vec3,PtNo,ta)
'  XShade.Message"x="&round(Pos(0)*1000)/1000&_
'  "   y="&round(Pos(1)*1000)/1000&_
'  "   z="&round(Pos(2)*1000)/1000

'--------------------ベジェ曲線でta(始点)からtb(終点)までの長さを計算
'    (ｱﾝｶｰﾎﾟｲﾝﾄ番号PtNo〜PtNo+1の部分が対象）
'  n=20: ta=0: tb=1: PtNo=0
'  b=SimpLength(Vec3,n,ta,tb,PtNo)
'  XShade.Message"ta="&ta&"から tb="&tb&"までの距離は"&_
'                 (round(b*1000)/1000)&"です。"

'--------------------ベジェ曲線上の距離を指定してtﾊﾟﾗﾒｰﾀを計算する
'    (ｱﾝｶｰﾎﾟｲﾝﾄ番号PtNo〜PtNo+1の部分が対象）
'  PtNo=0: Ltarget=3000: n=10
'  tb=tParam(Vec3,PtNo,Ltarget,n)
'  if tb>=0 and tb<=1 then
'      XShade.Message"ｱﾝｶｰﾎﾟｲﾝﾄ番号="&PtNo&"から"&_
'      Ltarget&"mmの地点のtﾊﾟﾗﾒｰﾀは"&round(tb*10000)/10000&"です。"
'  end if

'-------------------ｘｙｚいずれかを指定しベジェ曲線上でどのｔを取るか計算
'    (ｱﾝｶｰﾎﾟｲﾝﾄ番号PtNo〜PtNo+1の部分が対象）
'  PtNo=0: DatXYZ=300: Direction=0
'  tt=XyzToT(PtNo,Vec3,DatXYZ,Direction)
'  if tt(0)<>false then
'    nd=ubound(tt,1)
'    XShade.Message"解の数は"&nd&"個です。"
'    for i=1 to nd
'      XShade.Message"解("&i&")は"&round(tt(i)*100000)/100000
'    next
'  end if


'==========================曲率関連=================================
'-------------------ベジェ曲線でｔを指定し、その点の曲率、曲率半径を計算
'    (ｱﾝｶｰﾎﾟｲﾝﾄ番号PtNo〜PtNo+1の部分が対象）
'  PtNo=0: tb=0.55
'  a=Kyokuritu(Vec3,PtNo,tb)

'-------------------ベジェ曲線でｔを指定し、その点の接線の単位ベクトルを計算
'    (ｱﾝｶｰﾎﾟｲﾝﾄ番号PtNo〜PtNo+1の部分が対象）
'  PtNo=0: tb=0.55
'  eVec=eVecSessen(Vec3,PtNo,tb)
'  if not IsArray(eVec) then exit sub

'-------------------ベジェ曲線でｔを指定し、その点の法線の単位ベクトルを計算
'    (ｱﾝｶｰﾎﾟｲﾝﾄ番号PtNo〜PtNo+1の部分が対象）
'  PtNo=0: tb=0.55
'  eVec=eVecHousen(Vec3,PtNo,tb)
'  if not IsArray(eVec) then exit sub

'-------------------ベジェ曲線でｔを指定し、その点で次法線の単位ベクトルを計算
'    (ｱﾝｶｰﾎﾟｲﾝﾄ番号PtNo〜PtNo+1の部分が対象）
'  PtNo=0: tb=0.55
'  eVec=eVecJihou(Vec3,PtNo,tb)
'  if not IsArray(eVec) then exit sub


'---単位ベクトルを計算したものの線分をShade上に描画する。
'  l=400
'  Pos=calcPos(Vec3,PtNo,tb)
'  Call testDrawLine(eVec,Pos,PtNo,tb,L)
End sub


'===================================================================
'=                  主要なプロシージャー(subまたはfunction)
'===================================================================

'-----------tパラメータから座標値（のみ）の計算
Function calcPos(Vec3,PtNo,ta)
'階差とtパラメータ(=ta)から座標値のみを求める。
'戻値---tパラメータtaでの座標値
'引数---
' Vec3  ：(in)座標ﾃﾞｰﾀ(ｱﾝｶｰﾎﾟｲﾝﾄ番号、座標軸方向[x,y,z]、ﾎﾟｲﾝﾄの種類[ｱﾝｶｰ,In,Out])
' PtNo  ：(in)取り出す座標値のｱﾝｶｰﾎﾟｲﾝﾄ番号(指定する)
' ta     ：(in)求めたいtパラメータ

  Dim XYZ0(2),XYZ1(2),XYZ2(2),XYZ3(2)  '階差
  Dim Pos(2)  '：座標値用配列
  Dim P       '：３次のベジエ制御点P0〜P3の座標値
  Dim j       '：カウンタ用変数(xyz)

  P=GetnVec(Vec3,PtNo)
  Call kaisa(P,XYZ0,XYZ1,XYZ2,XYZ3)

  For j = 0 To 2
    Pos(j) = XYZ3(j)*ta^3+3*XYZ2(j)*ta^2+3*XYZ1(j)*ta+XYZ0(j)
  Next
  calcPos=Pos
End function


'--------------ベジェ曲線の弧長（S）を求める
'ベジェパラメーターtaからtb間での長さを求めるFunction
Function SimpLength(Vec3,n,ta,tb,PtNo)
'引数---
'  Vec3：(in)線型状の座標ﾃﾞｰﾀ【2次元】が格納されている。
'  n   ：(in)分割数(精度に関係する)
'  ta  ：(in)最初のﾎﾟｲﾝﾄのベジェパラメータ
'  tb  ：(in)終点のベジェパラメータ
'  PtNo：(in)PtNo番目のｱﾝｶｰﾎﾟｲﾝﾄ〜PtNo+1番目のｱﾝｶｰポイントまで計算される
  dim XYZ0(3),XYZ1(3),XYZ2(3),XYZ3(3)
  Dim Pos(3), dPos(3), ddPos(3)
  Dim P
  Dim j,k,Ft1,Ft2n,Ft,SFt1,SFt2
  Dim A,B,C,h,t()
  '---引数ｎが決まったためｔを動的に宣言
  Redim t(2*n)

  'ベジェ計算用の座標ﾃﾞｰﾀを取り出す
  P=GetnVec(Vec3,PtNo)
  '---ベジェパラメータｔを2n等分にしそれぞれのｋ番目のベジェパラメータを配列に割り付け
  for k=0 to 2*n
    t(k)=ta+k*(tb-ta)/(2*n)
  next

  '-----Simpson法で必要な各パラメータの算出-----
  Call kaisa(P,XYZ0,XYZ1,XYZ2,XYZ3)                     '階差の計算
  Call tBibun(XYZ0,XYZ1,XYZ2,XYZ3,t(0),Pos,dPos,ddPos)  'XYZ,1階微分,2階微分を計算
  Ft1=SQR(dPos(0)^2+dPos(1)^2+dPos(2)^2)
  Call kaisa(P,XYZ0,XYZ1,XYZ2,XYZ3)
  Call tBibun(XYZ0,XYZ1,XYZ2,XYZ3,t(2*n),Pos,dPos,ddPos)
  Ft2n=SQR(dPos(0)^2+dPos(1)^2+dPos(2)^2)

  SFt1=0
  for k=1 to n
    Call kaisa(P,XYZ0,XYZ1,XYZ2,XYZ3)
    Call tBibun(XYZ0,XYZ1,XYZ2,XYZ3,t(2*k-1),Pos,dPos,ddPos)
    Ft=SQR(dPos(0)^2+dPos(1)^2+dPos(2)^2)
    SFt1=SFt1+Ft
  next

  SFt2=0
  for k=1 to n-1
    Call kaisa(P,XYZ0,XYZ1,XYZ2,XYZ3)
    Call tBibun(XYZ0,XYZ1,XYZ2,XYZ3,t(2*k),Pos,dPos,ddPos)
    Ft=SQR(dPos(0)^2+dPos(1)^2+dPos(2)^2)
    SFt2=SFt2+Ft
  next

  'シンプソンの公式に必要なパラメータ
  A=Ft1+Ft2n	        'パラメータA
  h=(tb-ta)/(2*n)	'パラメータh
  B=4*SFt1		'パラメータB
  C=2*SFt2		'パラメータC

  'Simpsonの公式
  SimpLength=(h/3)*(A+B+C)
  XShade.Message PtNo&"番目のポイントと"&PtNo+1&"番目のポイントの長さ：S="&Round(SimpLength)
End function


'--------------tパラメータ算出（指定ﾎﾟｲﾝﾄからの距離でtパラメーターを計算）
Function tParam(Vec3,PtNo,Ltarget,n)
'ｱﾝｶｰﾎﾟｲﾝﾄからLtargetの距離にある曲線上の点のtパラメータを計算する
'ついでにそのﾎﾟｲﾝﾄのｘｙｚ座標も参考までに計算
'戻値---tパラメーター
'引数---
' Vec3   ：(in)対象となる線形状の座標データ(3次元）
' PtNo   ：(in)PtNo番目のｱﾝｶｰﾎﾟｲﾝﾄ〜PtNo+1番目のｱﾝｶｰポイントまでが計算の対象となる。
' Ltarget：(in)アンカーポイントからの距離（指定値)
' n      ：(in)分割数 弧長計算(SimpLength)で使用する。（精度に影響する）

  XShade.Message"Ltarget="&Ltarget
  Dim XYZ0(2),XYZ1(2),XYZ2(2),XYZ3(2)  '階差
  Dim ta,tb          'ta：開始tパラメータ  tb：終了tパラメータ
  Dim dLnew          'dLnew：計算値と目標値との差(長さ)
  Dim flag,alpha     'flag：判定用フラグ  alpha：収束判定値(許容誤差範囲)
  Dim P               '３次のベジエ制御点P0〜P3の座標値
  Dim Pos(2)         'ベジェ曲線上の点の座標値
  Dim dPos(2)        'Posの微分
  Dim ddPos(2)       'Posの二階微分

  ta=0:tb=1:alpha=0.0001

  P=GetnVec(Vec3,PtNo)   '３次のベジエ制御点P0〜P3の座標値を取得する
  tb=LTarget/SimpLength(Vec3,n,ta,tb,PtNo)   '弧長を求める		'
  if abs((SimpLength(Vec3,n,ta,tb,PtNo)-Ltarget))/Ltarget&gtalpha then
    Flag=0
    '---収束条件に当てはまるまで逐次計算
    Do
      Call kaisa(P,XYZ0,XYZ1,XYZ2,XYZ3)   '階差を求める	
      Call tBibun(XYZ0,XYZ1,XYZ2,XYZ3,tb,Pos,dPos,ddPos) 'Pos,Posの微分,二階微分を求める
      dLnew=SQR(dPos(0)^2+dPos(1)^2+dPos(2)^2)
      tb=tb+(Ltarget-SimpLength(Vec3,n,ta,tb,PtNo))/dLnew
      if abs(SimpLength(Vec3,n,ta,tb,PtNo)-LTarget)/LTarget < alpha then flag=1
    Loop until Flag
  end if
  tParam=tb
  Call tBibun(XYZ0,XYZ1,XYZ2,XYZ3,tb,Pos,dPos,ddPos)
End function


'-----------XYZのうち一つの値を渡し、ベジェ曲線と交差する交点のtﾊﾟﾗﾒｰﾀを計算する(誤差あり）
Function XyzToT(PtNo,Vec3,Dat,j)
'戻値---tパラメータ
'引数---
' PtNo：(in)PtNo番目のｱﾝｶｰﾎﾟｲﾝﾄ〜PtNo+1番目のｱﾝｶｰポイントまでが計算の対象となる。
' Vec3：(in)対象となる線形状の座標データ(3次元配列）
' Dat ：(in)次のjで指定した座標軸方向の数値(XYZのうち一つの値)
' j   ：(in)指定する座標軸方向（0：x , 1：y , 2：z）
  Dim XYZ0(2),XYZ1(2),XYZ2(2),XYZ3(2)  '階差
  Dim P                  '３次のベジエ制御点P0〜P3の座標値
  Dim a,b,c              'カルダノ法用パラメータ(3次式の項の係数)
  Dim tt                 'tパラメータ(データは複数(1〜3)の場合もある）
  Dim tt2()              '範囲内のtパラメータ
  Dim nd                 '返されるtパラメータの数
  Dim nd2                '範囲内の解(tパラメータ)の数
  Dim i

  nd2=0
  redim tt2(nd2)

  P=GetnVec(Vec3,PtNo)  'ベジエ制御点P0〜P3の座標値を取り出す

  Call kaisa(P,XYZ0,XYZ1,XYZ2,XYZ3)  '1〜3段階差を求る

  if XYZ3(j)<>0 then           '各次元の係数で処理を振り分ける。
    a=3*XYZ2(j)/XYZ3(j)
    b=3*XYZ1(j)/XYZ3(j)
    c=(XYZ0(j)-Dat)/XYZ3(j)
    tt=Cardano(a,b,c)          '３次式の解を求める
  elseif XYZ2(j)<>0 then
   a=3*XYZ2(j)
    b=3*XYZ1(j)
    c=(XYZ0(j)-Dat)
    tt=Formula2D(a,b,c)        '２次式の解を求める
  elseif XYZ1(j)<>0 then
    a=3*XYZ1(j)
    b=(XYZ0(j)-Dat)
    tt=Formula1D(a,b)          '１次式の解を求める
  else
    XShade.message"この方向には変化量がありません"
    tt2(0)=false
    exit function
  end if

  nd=Ubound(tt,1)
  if tt(0)=false then
    XShade.message"交差ポイントはありません"
    tt2(0)=false
    exit function
  end if

  '範囲(0〜1)外の解を排除
  tt2(0)=true
  for i=1 to nd
    if tt(i)>=0 and tt(i)<=1 then 
      nd2=nd2+1
      redim preserve tt2(nd2)
      tt2(nd2)=tt(i)
    end if
  next
  if nd2=0 then         '解が０の場合はfalseを返す。
    tt2(0)=false
    exit function
  end if

  XyzToT=tt2
end Function


'--------------ベジェ曲線の曲率（κ）を求める
Function kyokuritu(Vec3,PtNo,tb)
'戻値---曲率κ
'引数---
' Vec3：(in)対象となる線形状の座標データ(3次元）
' PtNo：(in)PtNo番目のｱﾝｶｰﾎﾟｲﾝﾄ〜PtNo+1番目のｱﾝｶｰポイントまでが計算の対象となる。
' tb  ：(in)tパラメータ(指定値）
' Pos ：(out)曲線上のtb点(曲率を求めた点)の座標
  Dim XYZ0(2),XYZ1(2),XYZ2(2),XYZ3(2)  '階差
  Dim Kapper,Rou          'Kapper：曲率  Rou：曲率半径
  Dim Pos(2), dPos(2), ddPos(2)  'Posの微分、Posの二階微分(tで微分)
  Dim P                       '３次のベジエ制御点P0〜P3の座標値
  Dim sPos(2),ssPos(2)   'Posの微分、Posの二階微分(s弧長で微分)

  'ベジェ計算用の座標ﾃﾞｰﾀを取り出す
  P=GetnVec(Vec3,PtNo)
  '階差を求めｔの一階微分、二階微分を計算する
  Call kaisa(P,XYZ0,XYZ1,XYZ2,XYZ3)
  Call tBibun(XYZ0,XYZ1,XYZ2,XYZ3,tb,Pos,dPos,ddPos)
  'ｓの一階微分、二階微分を計算する
  Call sBibun(dPos,ddPos,sPos,ssPos)
  'κ(曲率）を求める公式
  kyokuritu=sqr(ssPos(0)^2+ssPos(1)^2+ssPos(2)^2)
  Rou=1/kyokuritu
  XShade.Message PtNo&"番目で、ベジェパラメータt="&tb&"のときの曲率κ="_
            &round(kyokuritu*1000000)/1000000
  XShade.Message PtNo&"番目で、ベジェパラメータt="&tb&"のときの曲率半径ρ="_
            &round(Rou*1000)/1000
End Function


'--------------ベジェ曲線の接線パラメーター
Function eVecSessen(Vec3,PtNo,tb)
'ベジェパラメータｔでのベジェ曲線の接線のﾊﾟﾗﾒｰﾀeVec(l,m,n)を算出する

'戻値---接線の方向ベクトル
'引数---
' PtNo：(in)PtNo番目のｱﾝｶｰﾎﾟｲﾝﾄ〜PtNo+1番目のｱﾝｶｰポイントまでが計算の対象となる。
' tb  ：(in)tパラメータ(指定値）
' Vec3：(in)対象となる線形状の座標データ(3次元配列）
' Pos ：(out)曲線上のtb点(接線を求めた点)の座標
  Dim XYZ0(2),XYZ1(2),XYZ2(2),XYZ3(2)  '階差
  Dim Kapper,Rou         'Kapper：曲率  Rou：曲率半径
  Dim Pos(2), dPos(2), ddPos(2)  'Posの微分、Posの二階微分(tで微分)
  Dim P                       '３次のベジエ制御点P0〜P3の座標値
  Dim sPos(2),ssPos(2)   'Posの微分、Posの二階微分(s弧長で微分)
  Dim eVec(2)            '接線の方向(単位)ベクトル
  Dim j,S

  P=GetnVec(Vec3,PtNo)
  '階差を求めｔの一階微分、二階微分を計算する
  Call kaisa(P,XYZ0,XYZ1,XYZ2,XYZ3)
  Call tBibun(XYZ0,XYZ1,XYZ2,XYZ3,tb,Pos,dPos,ddPos)
  S=sqr(dPos(0)^2+dPos(1)^2+dPos(2)^2)
  If S=0 then
    eVecSessen=false:Exit Function
  end if
  for j=0 to 2
    eVec(j)=dPos(j)/S
  next
  eVecSessen=eVec
End Function


'--------------ベジェ曲線の法線パラメーター
Function eVecHousen(Vec3,PtNo,tb)
'ベジェパラメータｔでのベジェ曲線の法線のパラメータeVec(l,m,n)を算出する

'戻値---法線の方向ベクトル
'引数---
' PtNo：(in)PtNo番目のｱﾝｶｰﾎﾟｲﾝﾄ〜PtNo+1番目のｱﾝｶｰポイントまでが計算の対象となる。
' tb  ：(in)tパラメータ(指定値）
' Vec3：(in)対象となる線形状の座標データ(3次元配列）
' Pos ：(out)曲線上のtb点(法線を求めた点)の座標
  Dim XYZ0(2),XYZ1(2),XYZ2(2),XYZ3(2)  '階差
  Dim Kapper,Rou          'Kapper：曲率  Rou：曲率半径
  Dim Pos(2), dPos(2), ddPos(2)  'Posの微分、Posの二階微分(tで微分)
  Dim P                       '３次のベジエ制御点P0〜P3の座標値
  Dim sPos(2),ssPos(2)   'Posの微分、Posの二階微分(s弧長で微分)
  Dim eVec(2)              '法線の方向(単位)ベクトル
  Dim j

  P=GetnVec(Vec3,PtNo)
  '階差を求めｔの一階微分、二階微分を計算する
  Call kaisa(P,XYZ0,XYZ1,XYZ2,XYZ3)
  Call tBibun(XYZ0,XYZ1,XYZ2,XYZ3,tb,Pos,dPos,ddPos)
  'ｓの一階微分、二階微分を計算する
  Call sBibun(dPos,ddPos,sPos,ssPos)
  'κ(曲率）を求める公式
  Kapper=sqr(ssPos(0)^2+ssPos(1)^2+ssPos(2)^2)
  If Kapper=0 then
    eVecHousen=false:Exit Function
  end if
  if Kapper then
    Rou=1/Kapper
    for j=0 to 2
      eVec(j)=Rou*ssPos(j)	
    next
  end if
  eVecHousen=eVec
End Function 


'--------------ベジェ曲線の次法線パラメーター
Function eVecJihou(Vec3,PtNo,tb)
'ベジェパラメータｔでのベジェ曲線の法線のパラメータ(l,m,n)を算出する

'戻値---法線の方向ベクトル
'引数---
' PtNo：(in)PtNo番目のｱﾝｶｰﾎﾟｲﾝﾄ〜PtNo+1番目のｱﾝｶｰポイントまでが計算の対象となる。
' tb  ：(in)tパラメータ(指定値）
' Vec3：(in)対象となる線形状の座標データ(3次元配列）
' Pos ：(out)曲線上のtb点(法線を求めた点)の座標
  Dim XYZ0(2),XYZ1(2),XYZ2(2),XYZ3(2)  '階差
  Dim Kapper,Rou          'Kapper：曲率  Rou：曲率半径
  Dim Pos(2), dPos(2), ddPos(2)  'Posの微分、Posの二階微分(tで微分)
  Dim P                       '３次のベジエ制御点P0〜P3の座標値
  Dim sPos(2),ssPos(2)   'Posの微分、Posの二階微分(s弧長で微分)
  Dim eVec(2)              '法線の方向(単位)ベクトル
  Dim j
  P=GetnVec(Vec3,PtNo)
    '階差を求めｔの一階微分、二階微分を計算する
  Call kaisa(P,XYZ0,XYZ1,XYZ2,XYZ3)
  Call tBibun(XYZ0,XYZ1,XYZ2,XYZ3,tb,Pos,dPos,ddPos)
  'ｓの一階微分、二階微分を計算する
  Call sBibun(dPos,ddPos,sPos,ssPos)
  'κ(曲率）を求める公式
  Kapper=sqr(ssPos(0)^2+ssPos(1)^2+ssPos(2)^2)
  If Kapper=0 then
    eVecJihou=false:Exit Function
  end if
  if Kapper then
    Rou=1/Kapper
    eVec(0)=Rou*(sPos(1)*ssPos(2)-sPos(2)*ssPos(1))	
    eVec(1)=Rou*(sPos(2)*ssPos(0)-sPos(0)*ssPos(2))
    eVec(2)=Rou*(sPos(0)*ssPos(1)-sPos(1)*ssPos(0))		
  end if
  eVecJihou=eVec
End Function


'----------------------------描画ﾃｽﾄプロシージャ
sub testDrawLine(eVec,Pos,PtNo,tb,l)
'計算された接線、法線、次法線を確認するためにtest描画する
'戻値---無し
'引数---
' eVec   ：(in)方向ベクトル
' Pos    ：(in)描画開始ポイントの座標値
' PtNo   ：(in)ｱﾝｶｰﾎﾟｲﾝﾄの番号
' tb     ：(in)ベジェパラメータ(値を表示するだけに使う)
' l      ：(in)描画する線分の長さ

  '除算で分母がゼロになる場合の処理
  if eVec(0) or eVec(1) or eVec(2) then
    XShade.Message"eVec=( "&round(eVec(0)*100)/100&" , "_
     &round(eVec(1)*100)/100&" , "&round(eVec(2)*100)/100&" )"
    XShade.Message"Pos=( "&round(Pos(0)*100)/100&" , "_
     &round(Pos(1)*100)/100&" , "&round(Pos(2)*100)/100&" )"
  else
    XShade.Message"方向ベクトルがゼロベクトルなので描画できません。"
  end if

  '曲線上の指定されたベジェパラメータ位置にﾎﾟｲﾝﾄを追加する
  XShade.InsertControlPoint(PtNo+tb)

  'ラインを描画する
  XShade.StartPolygon(0)
    XShade.AppendPoint Pos(0),Pos(1),Pos(2)
    XShade.AppendPoint Pos(0)+l*eVec(0),Pos(1)+l*eVec(1),Pos(2)+l*eVec(2)
  XShade.FinishPolygon()
  XShade.SelectSister 1
End Sub

'================================================================================
'=     子プロシージャ-の記述　
'================================================================================

'-----------座標値,微分,二階微分の計算
Sub tBibun(XYZ0,XYZ1,XYZ2,XYZ3,ta,Pos,dPos,ddPos)
'階差とtパラメータ(=ta)から座標値,微分,二階微分を求める。
'戻値---無し
'引数---
' XYZ0   ：(in)P0座標値(1次元、3要素の配列）
' XYZ1   ：(in)第一段階差(1次元、3要素の配列）
' XYZ2   ：(in)第二段階差(1次元、3要素の配列）
' XYZ3   ：(in)第三段階差(1次元、3要素の配列）
' ta     ：(in)求めたいtパラメータ
' Pos    ：(out)tでのxyz座標値
' dPos   ：(out)tでの微分値
' ddPos  ：(out)tでの二階微分値

  Dim j   '：カウンタ用変数(xyz)

  'XYZ座標､一階微分､二階微分
  For j = 0 To 2
    Pos(j) = XYZ3(j)*ta^3+3*XYZ2(j)*ta^2+3*XYZ1(j)*ta+XYZ0(j)
    dPos(j) = 3*(XYZ3(j)*ta^2+2*XYZ2(j)*ta+XYZ1(j))
    ddPos(j) = 6*(XYZ3(j)*ta+XYZ2(j))
  Next
End Sub

'-------------ｔの微分，t2の階微分から弧長の微分、弧長の2階微分を計算
Function sBibun(d1,d2,S,SS)
'戻値---成功か失敗か(成功：true、失敗：false）
'引数---
' d1  ：(in)t微分
' d2  ：(in)ｔ二階微分
' S   ：(out)s微分
' SS  ：(out)s二階微分
  Dim A,B,j
  A=d1(0)^2+d1(1)^2+d1(2)^2
  B=d1(0)*d2(0)+d1(1)*d2(1)+d1(2)*d2(2)
  if A=0 or B=0 then
    sBibun=false
    exit function
  end if

  for j=0 to 2
    S(j)=d1(j)/sqr(A)
    SS(j)=(d2(j)/A)-(d1(j)*B/(A^2))
  next
  sBibun=true
End Function

'---選択している線型状の座標ﾃﾞｰﾀを取得する(３次元)
Function GetVec3(NA)
'戻値---ｱﾝｶｰﾎﾟｲﾝﾄの座標ﾃﾞｰﾀ　3次元配列
'　　　　各次元は (ｱﾝｶｰﾎﾟｲﾝﾄ番号、座標軸方向[x,y,z]、ﾎﾟｲﾝﾄの種類[ｱﾝｶｰ,In,Out])
'引数---
' NA   ：アンカーポイントの数
' 
  Dim Vec()          'ｱﾝｶｰﾎﾟｲﾝﾄの座標ﾃﾞｰﾀ(データ構造は戻り値と同じ)
  Dim i,j            'for next用カウンタ
  Redim Vec(NA,3,3)
  for i=0 to NA-1
    for j=0 to 2
      Vec(i,j,0)=XShade.AnchorPoint(i,j)
      Vec(i,j,1)=XShade.InHandle(i,j)
      Vec(i,j,2)=XShade.OutHandle(i,j)
    next
  next
  GetVec3=Vec
End Function

'------------------------------ベジェ計算に使うP0〜P3の座標をVecから拾い出す	
Function GetnVec(Vec3,PtNo)
'戻値---'３次のベジエ制御点P0〜P3の座標値(2次元)
'引数---
' Vec3  ：(in)座標ﾃﾞｰﾀ(ｱﾝｶｰﾎﾟｲﾝﾄ番号、座標軸方向[x,y,z]、ﾎﾟｲﾝﾄの種類[ｱﾝｶｰ,In,Out])
' PtNo  ：(in)取り出す座標値のｱﾝｶｰﾎﾟｲﾝﾄ番号(指定する)
  Dim j,P(3,2)
  for j=0 to 2 
    P(0,j)=Vec3(PtNo,j,0)                'P0:最初のポイントの座標値		
    P(1,j)=Vec3(PtNo,j,2)                'P1:最初のポイントのInハンドル座標値
    P(2,j)=Vec3(PtNo+1,j,1)              'P2:次のポイントのOutハンドル座標値
    P(3,j)=Vec3(PtNo+1,j,0)              'P3:次のポイントの座標値
  next
  GetnVec=P
End Function


'--------------１次式の解法
function Formula1D(a,b)
'１次式の解を求める。
'戻値---解(tパラメータ)　配列(1次元）0番目は成否に使用する。
'引数---
' a   ：(in)1次式の１次の項の係数。
' b   ：(in)1次式の定数

  Dim t(1)                '解（ｔパラメータ）用変数

  t(0)=true
  t(1)=-b/a
  Formula1D=t
end function

'--------------２次式の解法(解の公式)
function Formula2D(a,b,c)
'２次式の解を求める。
'戻値---解(tパラメータ)　配列(1次元）0番目は成否に使用する。
'引数---
' a   ：(in)２次式の２次の項の係数
' b   ：(in)２次式の１次の項の係数
' c   ：(in)２次式の定数項

  Dim D                  '判別式の結果用の変数
  Dim nd                 '解の数を保管する変数
  Dim t()                '解（ｔパラメータ）用変数

  D=round((b^2-4*a*c)*1000000)/1000000
  if D>0 then
    nd=2
    Redim t(nd)
    t(0)=true
    t(1)=(-b+sqr(b^2-4*a*c))/(2*a)
    t(2)=(-b-sqr(b^2-4*a*c))/(2*a)
  elseif D=0 then
    nd=1
    redim t(nd)
    t(0)=true
    t(1)=(-b+sqr(b^2-4*a*c))/(2*a)
  elseif D<0 then
    nd=0
    Redim t(nd)
    t(0)=false
  end if

  Formula2D=t
end function

'--------------３次式の解法（カルダノ方式）
Function Cardano(a,b,c)
'３次式の解をカルダノの公式で求める。
'戻値---解(tパラメータ)　配列(1次元）0番目は成否に使用する。
'引数---
' a   ：(in)３次式の２次の項の係数
' b   ：(in)３次式の１次の項の係数
' c   ：(in)３次式の定数項

  Dim D                  '判別式の結果用の変数
  Dim q,p,u,v,ux         '一時変数
  Dim nd                 '解の数を保管する変数
  Dim tt()               '解（ｔパラメータ）用変数

  p=-((a/3)^2)+(b/3)
  q=((a/3)^3)-(a/3)*(b/2)+(c/2)
  D=round(((q^2)+(p^3))*1000000)/1000000

  if D>0 then        '解が1つの場合
    nd=1
    Redim tt(nd)
    tt(0)=true
    if -q+sqr(q^2+p^3)>=0 then
      u=(-q+sqr(q^2+p^3))^(1/3)
    else
      u=-((-(-q+sqr(q^2+p^3)))^(1/3))
    end if
    if -q-sqr(q^2+p^3)>=0 then
      v=(-q-sqr(q^2+p^3))^(1/3)
    else
      v=-((-(-q-sqr(q^2+p^3)))^(1/3))
    end if
    tt(1)=u+v-a/3
  elseif D=0 then    '解が2つの場合
    nd=2
    Redim tt(nd)
    tt(0)=true
    if (-q)>=0 then
      tt(1)=2*(-q)^(1/3)-a/3
      tt(2)=-((-q)^(1/3))-a/3
    else
      tt(1)=-2*(-q)^(1/3)-a/3
      tt(2)=-(-q)^(1/3)-a/3
    end if
  elseif D<0 then    '解が3つの場合
    nd=3
    Redim tt(nd)
    tt(0)=true
    ux=q/(p*(-p)^(1/2))
    u=XShade.Acos(ux)
    tt(1)=2*sqr(-p)*cos(u/3)-a/3
    tt(2)=2*sqr(-p)*cos((u/3)+(2*3.14159265358979)/3)-a/3
    tt(3)=2*sqr(-p)*cos((u/3)+(4*3.14159265358979)/3)-a/3
  end if

  Cardano=tt
End Function

'-----------階差の計算（ベジエ制御点Pから階差を計算）
Sub kaisa(P,XYZ0,XYZ1,XYZ2,XYZ3)
'３次のベジエ制御点Pの座標値から１〜３段階差を計算する。
'戻値---無し
'引数---
' P      ：(in)ベジエ制御点Pの座標値(1次元、3要素の配列）
' XYZ0   ：(out)P0座標値(1次元、3要素の配列）
' XYZ1   ：(out)第一段階差(1次元、3要素の配列）
' XYZ2   ：(out)第二段階差(1次元、3要素の配列）
' XYZ3   ：(out)第三段階差(1次元、3要素の配列）

  Dim DVec11(2), DVec12(2), DVec13(2)   '：第一段階差用変数
  Dim DVec21(2), DVec22(2)              '：第二段階差用変数
  Dim DVec31(2)                         '：第三段階差用変数
  Dim j                                 '：カウンタ用変数(xyz)

  '第一段階差
  For j = 0 To 2
    DVec11(j) = P(1, j) - P(0, j)
    DVec12(j) = P(2, j) - P(1, j)
    DVec13(j) = P(3, j) - P(2, j)
  Next

  '第二段階差
  For j = 0 To 2
    DVec21(j) = DVec12(j) - DVec11(j)
    DVec22(j) = DVec13(j) - DVec12(j)
  Next

  '第三段階差
  For j = 0 To 2
    DVec31(j) = DVec22(j) - DVec21(j)
  Next

  '階差表示の定義
  For j = 0 To 2
    XYZ0(j) = P(0, j)
    XYZ1(j) = DVec11(j)
    XYZ2(j) = DVec21(j)
    XYZ3(j) = DVec31(j)
  Next
End sub

個々のサブルーチンの詳細を見てください。