# パラメトリック表示設定
set parametric

# 格子設定
set grid

# グラフ縦横比設定
set size ratio -1

# アルキメデス曲線の始点角θ0、終点角θ1の設定 [rad]
t0=8*pi   # 積分区間の下限
t1=100   # 積分区間の上限

# サンプル数設定
N=10000   # 積分区間の分割数

# 内端の回転角設定 a [rad]
a=5.5

# ヒゲゼンマイの全長 L
L=0.5*(t1*sqrt(1+t1**2)+log(t1+sqrt(1+t1**2))\
   -t0*sqrt(1+t0**2)-log(t0+sqrt(1+t0**2)))

# 被積分関数 
f(t)=h*(t*cos(t)*sin(Q)+t*sin(t)*cos(Q)) #実数部被積分関数
g(t)=-h*(t*cos(t)*cos(Q)-t*sin(t)*sin(Q)) #虚数部被積分関数

# 積分範囲の設定
set trange [t0:t1]
set samples N+1

# いったん別ファイルにデータを書き出す（変形後のx,y座標）
set table "table10.txt"

plot lastx=0.0,lastyr=0.0,lastyi=0.0,integr = 0.0,\
      integi = 0.0,old1=0,old2=0,old3=0,\integzr=0,\
      integzi=0,integzr2=0,integzi2=0,integL=0,zr=0,\
      zi=0,Qold=a,\
     "+" using 1 : (dx=$1-lastx, \
     s=0.5*($1*sqrt(1+$1**2)+log($1+sqrt(1+$1**2))\
        -t0*sqrt(1+t0**2)-log(t0+sqrt(1+t0**2))),\
     Q=a*(L-s)/L\+0.0001211*($1-t0)*\
         (sin(0.909341*$1+17.6221)-1.06),\
     h=(Q-Qold)/dx,\
     integr = ($0==0 ? 0.0:integr+dx*(f($1)+lastyr)*0.5), \
     integi = ($0==0 ? 0.0:integi+dx*(g($1)+lastyi)*0.5), \
     lastx=$1, \
     lastyr=f($1),\
     lastyi=g($1),\
     zr=integr+$1*cos($1)*cos(Q)-$1*sin($1)*sin(Q),\
     zi=integi+$1*cos($1)*sin(Q)+$1*sin($1)*cos(Q),\
     dzx=($0==0 ? 0.0 : zr-old1),\
     dzy=($0==0 ? 0.0 : zi-old2),\
     dL=sqrt(dzx**2+dzy**2),\
     integzr=integzr+(zr+old1)*0.5*dL,\
     integzi=integzi+(zi+old2)*0.5*dL,\
     integzr2=integzr2+(zr**2+old1**2)*0.5*dL,\
     integzi2=integzi2+(zi**2+old2**2)*0.5*dL,\
     integL=dL+old3,\
     old1=zr,\
     old2=zi,\
     old3=integL,\
     Qold=Q,\
     Ihx=integzi2/L-(integzi/L)**2,\
     Ihy=integzr2/L-(integzr/L)**2,\
     zr) : (zi) : (integL) : (integzr/L):(integzi/L) \
     w xyerrorbars

unset table

# グラフの線種設定
set style line 1 linetype 1
set style line 1 linecolor -1
set style line 1 linewidth 2

set style line 2 linetype 1
set style line 2 linecolor 1
set style line 2 linewidth 1

# x軸y軸の範囲設定
set xrange [-150:150]
set yrange [-150:150]

# x軸y軸のラベル設定
set xlabel "x/a"
set ylabel "y/a"

# グラフの描画
plot "table10.txt" u 2:3 w l linestyle 1 notitle

# 重心座標の描画
replot integzr/L,integzi/L w p title "G(x,y)"

# ヒゲ玉の描画
replot t0* cos(t),t0*sin(t) w l linestyle 2 notitle

# 出力フォーマットとオプションの指定
set terminal jpeg

# 出力ファイル名の指定
set output "koteitanhige5.jpeg"

# グラフの再描画
replot

# フォーマットと出力のリセット
set output
set terminal win
[EOF]