* string01.f ギター弦振動をシミュレーションする（一本だけ）
*
* １．弦を適当な質点に分割して、それぞれの点の動きを計算する。
* ２．初期のみ弦の任意の一点に、弦に対して直角に外力を与える。
* ３．各質点にかかる力は弦の張力のみを考慮する。
* ４．弦の両端は固定端とする。（←次の課題？）
* ５．弦の伸びによる張力の変化を無視する。
* ６．弦の曲がりによる力（剛性）を無視する。
* ７．摩擦などその他の力を無視する。
* 
      integer imax
      real px(0:8000),py(0:8000)
      real fx(0:8000),fy(0:8000)
      real vx(0:8000),vy(0:8000)
      real fx1,fx2,fy1,fy2
      real length,tension,mass,stress,body
      real hit_pos,hit_pow,hit_dur
      real dt,l,dx,dy,dx1,dy1,dx2,dy2,t
      real pi,py_max
      open(30,file="output3.dat")
      pi=4.0*atan(1.0)
***
* px,py          各質点の位置(m)
* fx,fy          各質点にかかる力(N)
* vx,vy          各質点の変位速度(m/s)
*
* imax=100       弦の分割数：100分割≒約6mm間隔
* dt=0.000000025 計算刻み時間(sec)：25ナノ秒＝50MHz
* length=0.648   弦長(m)：0.648m≒25.5インチ
* tension=1000.0 張力(N)：1000N≒102kgf（←大きすぎ。3〜8kgwが妥当？）
* g_mass=0.002   質量(kg)：0.002kg=2g（実際は0.5gぐらいか？）
* hit_pos=0.08   弾く位置(m)
* hit_pow=20.0   弾く力(N)：20N≒2kgf（←？）
* hit_dur=0.1    弾く最長時間(sec)（実際にはそれより前に弾き終わる）
* s_pos1=0.03    ピックアップ1の位置(m)
* s_pos2=0.15    ピックアップ2の位置(m)
*
      imax=100
      dt=0.000000025
      length=0.648
      tension=50.0
      g_mass=0.0005
      hit_pos=0.08
      hit_pow=1.0
      hit_dur=0.100
      s_pos1=0.03
      s_pos2=0.15
*      neck=0.4
*      body=3.6
*      stress=20000.0
*
      mass=g_mass/real(imax)
      i0=int(hit_pos/length*real(imax)+0.5)
      i1=int(s_pos1/length*real(imax)+0.5)
      i2=int(s_pos2/length*real(imax)+0.5)
*
      do 100 i=0,imax,1
        px(i)=real(i)/real(imax)*length
        py(i)=0.0
        fx(i)=0.0
        fy(i)=0.0
        vx(i)=0.0
        vy(i)=0.0
100   continue
      py_max=-9999.
*
      do 200 ii=0,40000000,1 
        t=real(ii)*dt
*        l=0.0
*        do 210 i=1,imax,1
*          dx=px(i)-px(i-1)
*          dy=py(i)-py(i-1)
*          l=l+sqrt(dx*dx+dy*dy)
*210     continue 
        do 220 i=1,imax-1,1
          if((t.lt.hit_dur).and.(i.eq.i0)) then
            if(py(i).gt.py_max) then
              py_max=py(i)
            else
              if((py_max.gt.0.001).and.(py(i).lt.py_max*0.8)) hit_dur=t
            end if        
          end if
          dx1=px(i+1)-px(i)
          dy1=py(i+1)-py(i)
          dx2=px(i-1)-px(i)
          dy2=py(i-1)-py(i)
          fx1=tension*dx1/sqrt(dx1*dx1+dy1*dy1)
          fy1=tension*dy1/sqrt(dx1*dx1+dy1*dy1)
          fx2=tension*dx2/sqrt(dx2*dx2+dy2*dy2)
          fy2=tension*dy2/sqrt(dx2*dx2+dy2*dy2)
          fx(i)=fx1+fx2
          fy(i)=fy1+fy2
          if((t.lt.hit_dur).and.(i.eq.i0)) fy(i)=fy(i)+hit_pow
          vx(i)=vx(i)+fx(i)/mass*dt
          vy(i)=vy(i)+fy(i)/mass*dt
220     continue
*
*        dx1=px(1)-px(0)
*        dy1=py(1)-py(0)
*        fy1=tension*dy1/sqrt(dx1*dx1+dy1*dy1)
*        fy2=-py(0)*stress
*        fy(0)=fy1+fy2
*        vy(0)=vy(0)+fy(0)/body*dt
*
*        dx1=px(imax-1)-px(imax)
*        dy1=py(imax-1)-py(imax)
*        fy1=tension*dy2/sqrt(dx2*dx2+dy2*dy2)
*        fy2=-py(imax)*stress
*        fy(imax)=fy1+fy2
*        vy(imax)=vy(imax)+fy(imax)/neck*dt
*
        do 230 i=0,imax,1
          px(i)=px(i)+vx(i)*dt
          py(i)=py(i)+vy(i)*dt
230     continue
*
        if(mod(ii,40).eq.0) then
          write(30,1030) t,i0,i1,i2,0,imax
     &          ,py(i0),py(i1),py(i2),py(0),py(imax)
1030      format(f12.6,5i5,8f12.8)
        end if
*
*        write(*,1000) ii,t,i0,px(i0),py(i0),fx(i0),fy(i0)
*     &                ,fy(0),fy(imax),py(0),py(imax)
*1000     format(i8,f12.7,i6,4f12.6,2f20.8,2f20.8) 
200   continue
      close(30)
      end