program heli
c     seismometer with lacoste type suspension
c  Eingabedaten:
c     xu,yu: lage der festen federklemme (cm)
c     xo,yo: dasselbe fuer die bewegliche klemme (cm)
c     fl: federlaenge ohne kraft (cm)
c     step: Suchschritt in (cm)
c     fk: federkonstante (kp/cm)
c     hin: Rckstellmoment der Gelenke
c
      implicit real*8 (a-h,o-z)
      dimension gg(3),hh(3),f(5,2,5),g(5,2,5)
      common ax,pi4,fl,hstep,hin,fk,torad,dnull
      character*24 text
      dnull=0.d0
      pi4=datan2(1.d0,1.d0)
      torad=pi4/45.d0
      hstep=1.d-2
      open(15,file='heli.in')
      open(16,file='heli.out')
   10 read (15,*,end=14) xu,yu,xo,yo,fl,step,fk,hin,ax,text
    1 format(8f5.2,a)
      if(fk.eq.0.) stop
    3 r=dsqrt(xo*xo+yo*yo)
      write(16,4) text,xu,yu,xo,yo,fl,step,fk,hin,ax
      write(6,4) text,xu,yu,xo,yo,fl,step,fk,hin,ax
    4 format(//a24," eingabedaten:"/ 
     -"      xu      yu      xo      yo      fl    step      fk     hin
     -      ax"
     -/9f8.3//" koordinatenursprung im gelenk; winkel in grad ccw von de
     -r x-achse"/" variiert wird xu und yu, also der gestellfeste einspa
     -nnpunkt der feder"/" einheiten: hertz, kilopond und cm"/
     -" fl=Federlaenge ohne Kraft; step=Suchschritt; fk=Federkonstante"/
     -" hin=Rueckstellmoment der Gelenke; ax=Richtung der empf. Achse")
      ax=ax*torad
      do 5 l=1,5
      xuv=xu+(l-3)*step
      do 5 n=1,5
      yuv=yu+(3-n)*step
      call freq(xuv,yuv,wu,xo,yo,wo,fq,tq,gg,hh,ter,xk,yk,d,r)
      f(l,1,n)=fq
      f(l,2,n)=tq
      g(l,1,n)=ter
    7 g(l,2,n)=hh(3)
    5 continue
      write(16,12)
      write(6,12)
   12 format(/ "  eigenfrequenz**2                     erste ableitung"
     & /)
      write(16,6) f
      write(6,6) f
    6 format(5f7.3,3x,5f7.3)
      write(16,13)
      write(6,13)
   13 format( /"  zweite ableitung                     gewicht"/)
      write(16,9) g
      write(6,9) g
    9 format(5f7.3,3x,5f7.3)
      goto 10
   14 stop
      end

      subroutine freq(xu,yu,wu,xo,yo,wo,fm,fdif,gg,hh,ter,xk,yk,d,r)

      implicit real*8 (a-h,o-z)
      dimension e(5),fkr(4),fko(3),fq(3),gg(3),hh(3)
      common ax,pi4,fl,hstep,hin,fk,torad,dnull
      zpiq = 8.d0*pi4
      zpiq = zpiq*zpiq
      ter=dnull
c  berechnen des elastischen potentials in 5 benachbarten lagen
      do 30 l=1,5
      h = hstep*(3-l)
      call kreis(xo,yo,b,h,xz,yz,bz,r)
      call such(xk,yk,d,xz-xu,yz-yu,bz,e(l))
      if(l.ne.3) goto 30
   30 continue
c  drehmoment und seismische masse
      dreh=(e(4)-e(2))/(2.*hstep)
      rmg=dreh/dcos(ax)
      rm=rmg/r
c  hinzunahme des schwerepotentials
      do 71 l=1,5
   71 e(l)=e(l)+rmg*dsin(ax+hstep*(3-l))
      fm=e(3)
      do 31 l=1,5
   31 e(l)=e(l)-fm
c  drehmoment als ableitung des potentials nach dem drehwinkel
      do 70 l=1,4
   70 fkr(l) =(e(l+1)-e(l))/hstep
c  rueckstellmoment als ableitung des drehmoments
      do 80 l=1,3
      fko(l) =(fkr(l+1)-fkr(l))/hstep+hin
c  quadrat der eigenfrequenz
   80 fq(l) = fko(l)/rm*981./r**2/zpiq
      fm=fq(2)
      fdif=0.5d0*(fq(3)-fq(1))/(hstep*r)
      ter=(fq(3)-fq(2)*2.+fq(1))/(hstep*r)**2
      hh(3)=rm
      return
   40 fm=dnull
      fdif=dnull
      do 50 j=1,3
      gg(j)=dnull
   50 hh(j)=dnull
      return
      end

      subroutine kreis(xo,yo,ba,defl,x,y,b,r)

      implicit real*8 (a-h,o-z)
      common ax,pi4,fl,hstep,hin,fk,torad,dnull
      del = datan2(yo,xo) +defl
      x=r*dcos(del)
      y=r*dsin(del)
      return
      end

      subroutine such(xk,yk,d,xz,yz,bz,e)
      implicit real*8 (a-h,o-z)
      common ax,pi4,fl,hstep,hin,fk,torad,dnull
      al=sqrt((xk-xz)**2+(yk-yz)**2)
      e=fk/2.*(al-fl)**2
      end