c     geometry of astatic leaf-spring suspensions   nov. 2000 ew
c
      program leaffind
c
c     Version 5, Maerz 2003, IGPP La Jolla
c
c     interne parameter:
c     acc: genauigkeitsforderung bei der iteration (cm)
c     xinc: inkrement beim berechnen der partiellen ableitungen (cm)
c     hstep: winkelinkrement bei der berechnung der pot. energie (rad)
c     maxit: abbruch der iteration nach maxit aufrufen von blf
c
c     achtung: es kommen drei verschiedene normierungen des drehmoments
c     vor. routine blf rechnet mit dimensionslosem drehmoment. bei der
c     inversion nach den startparametern ist dieses durch 5 dividiert.
c     das physikalische drehmoment in kp*cm ergibt sich durch multi-
c     plikation mit faktor*2*pi/180. derselbe faktor ist bei der
c     kraft anzuwenden.
c     analog kommt auch die kraft in zwei verschiedenen normierungen vor.
c
c     Bug fixed Dec 2007:
c     the sign of nlin1 was wrong, the freq routine has been corrected.
c     The problem is related to defining h=(3-l)*hstep, h going down.

      implicit double precision (a-h,o-z)
      dimension gg(3),hh(6),f(2,3,2),a(3),b(3),c(3)
      character*87 line
      common acc,xinc,ax,pi4,halb,ngg,gl,ramin,ramax,fl,itera,hstep,
     -nix,bu,maxit,fakt,dfak,fekru,torad,dnull,itotal,itrace
      common /spring/ xua,yua,xmin,xmax,ymin,ymax,xs(130),ys(130),last

      itrace=0
      ng=128
c     subdivision of spring. note that dimension of xs, ys must be ng+2.
      strmax=280.d0
c   1/strmax is the maximum permitted strain of the spring material
      dnull=0.d0
      halb=5.d-1
      pi4=datan2(halb,halb)
      torad=pi4/45.d0
      acc=1.d-10
      xinc=1.d-4
      maxit=3333
      nix=0
      hstep=1.d-2
      step=1.d-2
      ngg=ng-1
      iter=17
      itera=0

      open(7,file='leaf.in')
      open(8,file='leaf.out')
      write(6,100)
      write(8,100)
  100 format(' This program computes the mechanical characteristics'
     &/' of a leaf-spring suspension and searches for a spring'
     &/' geometry with a long and constant free period.'/
     &/' Files:    leaf.in ==> leafspring ==> leaf.out and spring.out'/)

      write(6,140)
      write(8,140)
  140 format(' Program option controlled by the first input parameter:'
     &    /' iopt = 0, all done'
     &    /'      = 1, compute properties of given spring geometry'
     &    /'      = 2, find astatic geometry with constant frequency'
     &    /'      = 3, find same with nonsymmetric circular spring'
     &    /'      = 4, find same with symmetric noncircular spring'
     &    /' iopt = 3 and 4: xo,yo,wo are determined from the desired'
     &    /'        spring geometry, not taken from the input file'
     &    /' iopt = 9, ignore this input line'/)

   10 read(7,'(a)',err=50,end=50) line
      read(line,'(i2)') iopt
      if(iopt.eq.9) goto 10
      if(iopt.eq.0) then
        write(6,150)
  150   format(' Option iopt=0 or blank in col.1 - stop')
        stop
      endif
       read (line,1,err=50,end=50) xu,yu,wu,xo,yo,wo,fl,ax,fekru,
     & breite,dicke,emodul,xki,yki,di,unf,unl
    1 format(2x,17f5.3)
      if(iopt.le.4) goto 70
   50 write(6,60)
   60 format(' Input error. The file leaf.in must contain the following'
     &/' 18 parameters in one line in format (i1,1x,17f5.3):'/)
      write(6,40)
      stop
   70 write(6,160) iopt
      write(8,160) iopt
  160 format(' Option in effect: iopt=',i1)
      if(iopt.eq.0) stop
      write(6,4)     xu,yu,wu,xo,yo,wo,fl,ax,fekru,breite,dicke,emodul,
     & xki,yki,di,unf,unl
      write(8,4)     xu,yu,wu,xo,yo,wo,fl,ax,fekru,breite,dicke,emodul,
     & xki,yki,di,unf,unl
    4 format(/' Parameters read from file leaf.in:'//
     &'       xu      yu      wu      xo      yo      wo      fl      ax
     &'/ 1x,8f8.3//
     &'       fk      wi      th              em     xki     yki      di
     &'/ 1x,3f8.3,f16.3,3f8.3//
     &'      unf     unl'
     & / 1x,2f8.3/)
      itotal=0
      write(6,40)
      write(8,40)
   40 format(' iopt:       program option, see above'/
     &' xu, yu, wu: coordinates and angle of one end of the spring'/
     &' xo, yo, wo: same for the other end of the spring. coordinates in
     &'/'             cm from the hinge. angles in degrees ccw from the'
     & /'             x axis. spring may be rotated around the hinge.'
     & /' fl, ax, fk: length of the spring, angle of the axis of sensiti
     &vity'/'             against the vertical, initial curvature of the 
     & spring'/'             (end-to-end, without load) in degrees.
     &'/' wi, th, em: width, thickness, elast. modulus of the spring
     &'/' xki,yki,di: estimated spring force and moment
     &'/' unf,unl:    eigenfrequency-squared and its derivative w.r.t.
     &'/'             position due to causes unaccounted for in this pro
     &gram'/)

      raref=dicke*strmax/2.d0
c     the elastic modulus is converted into the leaf-bending modulus
      emodul=emodul*16./15.
      fakt=emodul*breite*dicke**3/24.d0
      dfak=fakt*2.d0*torad
      gl=fl/ng
      last=0

      if(iopt.eq.3) then
        wor=wu*torad
        sehn=fl*dsin(wor)/wor
        xo=xu+sehn
        yo=yu
        wo=-wu
        xki=dnull
        yki=dnull
        di=2.d0*dfak*wu/fl
        write(6,80)
        write(8,80)
   80   format(' The spring is assumed to be a circular arc with the'
     & /' convex side upward. Input parameters xo, yo, and wo are ignore
     &d'/' and replaced by values calculated from the spring geometry.'
     &/' The same applies to xki, yki, and di.'/)
      endif
      if(iopt.eq.4) then
        wor=wu*torad
        xo=-xu
        yo=yu
        wo=-wu
        write(6,81)
        write(8,81)
   81   format(' The spring is assumed to be a symmetric arc with the'
     & /' convex side upward. Input parameters xo, yo, and wo are ignore
     &d'/' and replaced by values calculated from the spring geometry.'
     &/' The same applies to xki, yki, and di.'/)
      endif

      write(6,2) breite,dicke*10.d0,emodul/100.d0
      write(8,2) breite,dicke*10.d0,emodul/100.d0
    2 format(' spring: width ',f5.1,' cm, thickness ',f5.2,' mm, emodul 
     &',f7.0,' kp/mm**2'/)

      ax=ax*torad
      fekru=fekru*torad/ng
      xk=xki/dfak
      yk=yki/dfak
      d=di/dfak/5.d0

      if(iopt.eq.1) goto 110
      write(6,*) 'Iterative search:'
      write(8,*) 'Iterative search:'
  302 nix=0

   90 format(/' itera      fak        abalt        abneu         xk
     &         yk          d'/)

      do 5 l=1,2
      xuv=xu+(l-1.5)*step
      do 15 n=1,2
      yuv=yu+(n-1.5)*step
      if(iopt.eq.3) then
        xo=xuv+sehn
        yo=yuv
      endif
      if(iopt.eq.4) then
      xo=-xuv
      yo=yuv
      endif

      r=dsqrt(xo*xo+yo*yo)
      if(itrace.eq.1) write(6,90)  

      call freq(xuv,yuv,wu,xo,yo,wo,fq,tq,gg,hh,ter,xk,yk,d,r)
      fq=fq+unf
      tq=tq+unl

      f(l,1,n)=fq
      f(l,2,n)=tq
      f(l,3,n)=ter
   15 continue
    5 continue
c     write(8,12)
   12 format(/' frequency^2              nonlin1              nonlin2')
c     write(8,6) f
    6 format(2f6.3,9x,2f6.3,9x,2f6.3)
c     write(8,13)
   13 format(/' max. radius          min: radius              maxmass')
c     write(8,9) h
    9 format(2f6.2,9x,2f6.2,9x,2f6.3)
c     write(8,14)
   14 format(/'      xforce               yforce      torque at xu,yu')
c     write(8,19) g
   19 format(2f6.2,9x,2f6.2,9x,2f6.2)
      do 301 j=1,3
      c(j)=(f(1,j,1)+f(2,j,1)+f(1,j,2)+f(2,j,2))/4.d0
      a(j)=(f(2,j,1)+f(2,j,2)-f(1,j,1)-f(1,j,2))/2./step
      b(j)=(f(1,j,2)+f(2,j,2)-f(1,j,1)-f(2,j,1))/2./step
  301 continue
      det=a(1)*b(2)-a(2)*b(1)    
      write(6,303) xu,yu,c(1),c(2)
      write(8,303) xu,yu,c(1),c(2)
  303 format(' xu=',f7.3,'    yu=',f7.3,'    fr^2=',f10.6,'    nlin=',
     &f10.6)
      xcor=(c(1)*b(2)-b(1)*c(2))/det
      ycor=(a(1)*c(2)-a(2)*c(1))/det
      eps=dsqrt(1.d0/(1.d0+(xcor**2+ycor**2)/1.d-1))
      xu=xu-xcor*eps
      yu=yu-ycor*eps
      iter=iter-1
      if(c(1)**2+c(2)**2.gt.1d-9.and.iter.gt.0) then
      goto 302
      endif         
  110 write(6,30)
      write(8,30)
   30 format(/' Final spring geometry. Units are cm, kp, kg. Mass is 5 c
     &m from hinge.'/)
      if(iopt.eq.3) then
        xo=xu+sehn
        yo=yu
      endif
      if(iopt.eq.4) then
      xo=-xu
      yo=yu
      endif

      r=dsqrt(xo*xo+yo*yo)
      if(itrace.eq.1) write(6,90)
      nix=0
      call freq(xu,yu,wu,xo,yo,wo,fq,tq,gg,hh,ter,xk,yk,d,r)
      fq=fq+unf
      tq=tq+unl
      hh(3)=raref
      hh(6)=hh(5)*(hh(2)/raref)**3
      write(6,401) xu,yu,wu
      write(8,401) xu,yu,wu
  401 format('      xu=',f7.3,'         yu=',f7.3,'         wu=',f7.3)
      write(6,402) xo,yo,wo
      write(8,402) xo,yo,wo
  402 format('      xo=',f7.3,'         yo=',f7.3,'         wo=',f7.3)
      write(6,403) fq,tq,ter
      write(8,403) fq,tq,ter
  403 format('    fr^2=',f7.3,'      nlin1=',f7.3,'      nlin2=',f7.3)
      write(6,404) hh
      write(8,404) hh
  404 format('    rmax=',f7.3,'       rmin=',f7.3,'       rlim=',f7.3/
     &       '    mass=',f7.3,'  =?=  mass=',f7.3,'    masslim=',f7.3)
      write(6,405) gg
      write(8,405) gg
  405 format('    xfrc=',f7.3,'       yfrc=',f7.3,'       torq=',f7.3)
      write(6,406) strmax
      write(8,406) strmax
  406 format(/' Explanation of symbols: xu, yu, wu, xo, yo, wo as above.
     &'/' fr^2, nlin1, nlin2: squared frequency and its first two'/' der
     &ivatives with respect to the mass position (length unit = cm)'
     & /' rmax, rmin: max. and min. radius of curvature of the spring'
     & /' rlim: min. radius of curv. permitted by max. strain 1/',f4.0/
     &' mass: weight of the mass in kp, calculated in two different ways
     &'/' masslim: weight with a spring of same width but maximum thickn
     &ess'/' xfrc, yfrc, torq: force and torque in both clamps (kp and k
     &p*cm)')

      call final(xu,yu,wu,xo,yo,wo,fq,tq,gg,hh,ter,xk,yk,d,r)
      write(6,'(/a)') ' Extremal coordinates of the spring'
      write(8,'(/a)') ' Extremal coordinates of the spring'
      write(6,102) xmin+xua,xmax+xua,ymin+yua,ymax+yua
      write(8,102) xmin+xua,xmax+xua,ymin+yua,ymax+yua
  102 format("    xmin=",f6.3,"    xmax=",f6.3,
     &       "    ymin=",f6.3,"    ymax=",f6.3/
     &       " Shape of spring is stored in file spring.out")
      open(9,file='spring.out')
      write(9,104) (xua+xs(j),yua+ys(j),j=1,ng+2)
  104 format(2f10.3)
      close(9)

      write(6,407) itotal
      write(8,407) itotal
  407 format(/' Routine BLF was called ',i6,' times.'
     &/' Use xfrc, yfrc, torq as initial values xki, yki, di.'//
     &' ============================================================='/)
      goto 10
      end

      subroutine final(xu,yu,wu,xo,yo,wo,fm,fdif,gg,hh,ter,xk,yk,d,r)
      implicit double precision (a-h,o-z)
      dimension e(5),fkr(4),fko(3),fq(3),gg(3),hh(6)
      common acc,xinc,ax,pi4,halb,ngg,gl,ramin,ramax,fl,itera,hstep,
     -nix,bu,maxit,fakt,dfak,fekru,torad,dnull,itotal,itrace
      common /spring/ xua,yua,xmin,xmax,ymin,ymax,xs(130),ys(130),last
      xua=xu
      yua=yu
      last=1
      b=wo/10.d0
      bu=wu/10.d0
      zpiq = (8.*pi4)**2
      rm=5.d0      
      h = dnull
      call kreis(xo,yo,b,h,xz,yz,bz,r)
      call such(xk,yk,d,xz-xu,yz-yu,bz,edum)
      return
      end

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

      implicit double precision (a-h,o-z)
      dimension e(5),fkr(4),fko(3),fq(3),gg(3),hh(6)
      common acc,xinc,ax,pi4,halb,ngg,gl,ramin,ramax,fl,itera,hstep,
     -nix,bu,maxit,fakt,dfak,fekru,torad,dnull,itotal,itrace
      common /spring/ xua,yua,xmin,xmax,ymin,ymax,xs(130),ys(130),last
      xua=xu
      yua=yu

c     achtung... winkeleinheit fuer b ist 10 grad
c     achtung... wahres drehmoment ist  5.*d [*dfak]

      if(itrace.eq.1) write(6,15) xk,yk,d
   15 format(' freq: xk=',f7.3,'   yk=',f7.3,'    d=',f7.3)
      b=wo/10.d0
      bu=wu/10.d0
      zpiq = (8.*pi4)**2
      rm=5.d0
      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))
      e(l)=e(l)*fakt
      if(nix.eq.0) goto 40
      if(l.ne.3) goto 30
      hh(1)=ramin
      hh(2)=ramax
      hh(4)=(5.d0*d+yk*xu-xk*yu)*dfak/dcos(ax)/rm
      gg(1)=xk*dfak
      gg(2)=yk*dfak
      gg(3)=5.d0*d*dfak
   30 continue
c  drehmoment und seismische masse
      dreh=(e(4)-e(2))/(2.*hstep)
      rmg=dreh/dcos(ax)
      gm=rmg/rm
      hh(5)=gm
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
c  quadrat der eigenfrequenz
   80 fq(l) = fko(l)/gm*981./rm**2/zpiq
      fm=fq(2)
      fdif=0.5d0*(fq(1)-fq(3))/(hstep*rm)
      ter=(fq(3)-fq(2)*2.+fq(1))/(hstep*rm)**2
      return
   40 fm=dnull
      fdif=dnull
      do 50 j=1,3
      gg(j)=dnull
      hh(j+3)=dnull
   50 hh(j)=dnull
      return
      end

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

      implicit double precision (a-h,o-z)
      common acc,xinc,ax,pi4,halb,ngg,gl,ramin,ramax,fl,itera,hstep,
     -nix,bu,maxit,fakt,dfak,fekru,torad,dnull,itotal,itrace
      del = datan2(yo,xo) +defl
      x=r*dcos(del)
      y=r*dsin(del)
      b = ba+defl/torad/1.d1
      return
      end

      subroutine such(xk,yk,d,xz,yz,bz,e)

      implicit double precision (a-h,o-z)
      dimension a(3,3),dif(3),xkor(3)
      common acc,xinc,ax,pi4,halb,ngg,gl,ramin,ramax,fl,itera,hstep,
     -nix,bu,maxit,fakt,dfak,fekru,torad,dnull,itotal,itrace
      ab(a1,a2,a3)=a1*a1+a2*a2+a3*a3
      abq=acc*acc
      q=xinc
      itotal=itotal+itera
      itera=0
      fak=1.d00
    7 if(nix.eq.0) call blf(xk,yk,d,x,y,b,e)
    6 dif(1)=xz-x
      dif(2)=yz-y
      dif(3)=bz-b
      abneu=ab(dif(1),dif(2),dif(3))
      if(nix.gt.0) goto 5
    2 call blf(xk+q,yk,d,a(1,1),a(2,1),a(3,1),e)
      call blf(xk,yk+q,d,a(1,2),a(2,2),a(3,2),e)
      call blf(xk,yk,d+q,a(1,3),a(2,3),a(3,3),e)
      a(1,1)=(a(1,1)-x)/q
      a(2,1)=(a(2,1)-y)/q
      a(3,1)=(a(3,1)-b)/q
      a(1,2)=(a(1,2)-x)/q
      a(2,2)=(a(2,2)-y)/q
      a(3,2)=(a(3,2)-b)/q
      a(1,3)=(a(1,3)-x)/q
      a(2,3)=(a(2,3)-y)/q
      a(3,3)=(a(3,3)-b)/q
    5 xkalt=xk
      ykalt=yk
      dalt=d
      abalt=abneu
      xalt=x
      yalt=y
      balt=b
    3 call matin(a,dif,abneu*fak,xkor)
      xk=xkalt+xkor(1)
      yk=ykalt+xkor(2)
      d= dalt+xkor(3)
      nix=1
      call blf(xk,yk,d,x,y,b,e)
      dif(1) = xz-x
      dif(2) = yz-y
      dif(3) = bz-b
      abneu=ab(dif(1),dif(2),dif(3))
      if(itrace.eq.1) write(6,4) itera,fak,abalt,abneu,xk,yk,d
    4 format(1x,i5,3x,f6.0,2(3x,d10.3),3(3x,f8.2))
      if(itera.lt.maxit) goto 11
      nix=0
      return
   11 if(abneu.le.abq) then
        call blf(xk,yk,d,x,y,b,e)
        return
      endif
      if(abneu.lt.abalt) goto 2
      if(itera.gt.2) goto 14
      nix=0
      goto 7
   14 fak=fak*32.d00
      dif(1)=xz-xalt
      dif(2)=yz-yalt
      dif(3)=bz-balt
      goto 3
      end

      subroutine matin(a,dif,wq,xkor)

      implicit double precision (a-h,o-z)
      dimension a(3,3),dif(3),xkor(3),b(3,3),c(3)
      dnull=0.d0
      do 1 j=1,3
      c(j)=dnull
      do 1 k=1,3
      c(j)=c(j)+a(k,j)*dif(k)
      b(j,k)=dnull
      if(j.eq.k) b(j,k)=wq
      do 1 i=1,3
    1 b(j,k)=b(j,k)+a(i,j)*a(i,k)
      call gaus3(b,c,xkor)
      return
      end

      subroutine gaus3(aik,rs,f)

      implicit double precision (a-h,o-z)
      dimension aik(3,3),rs(3),f(3),h(4),imax(3)
      dnull=0.d0
      do 1401 j=1,3
      aikmax=dnull
      do 1402 k=1,3
      h(k)=aik(j,k)
      if(dabs(h(k)).le.aikmax) go to 1402
      aikmax=dabs(h(k))
      index=k
 1402 continue
      h(4)=rs(j)
      do 1403 k=1,3
      q=aik(k,index)/h(index)
      do 1404 l=1,3
 1404 aik(k,l)=aik(k,l)-q*h(l)
 1403 rs(k)=rs(k)-q*h(4)
      do 1405 k=1,3
 1405 aik(j,k)=h(k)
      rs(j)=h(4)
 1401 imax(j)=index
      do 1406 j=1,3
      index=imax(j)
 1406 f(index)=rs(j)/aik(j,index)
      return
      end


      subroutine blf(xkk,ykk,dk,xx,yy,bb,ee)
      implicit double precision (a-h,o-z)
      logical fin
      common acc,xinc,ax,pi4,halb,ngg,gl,ramin,ramax,fl,itera,hstep,
     -nix,bu,maxit,fakt,dfak,fekru,torad,dnull,itotal,itrace
      common /spring/ xua,yua,xmin,xmax,ymin,ymax,xs(130),ys(130),last
      fin=(last.eq.1)
      itera=itera+1
      xk = xkk
      yk = ykk
      d = dk*5.d0
      w=torad*gl
      e=dnull
      b=10.d0*bu*w/gl
      x=halb*dcos(b)*gl
      y=halb*dsin(b)*gl
      rmin=1.d12
      rmax=dnull
      if(fin) then
        xs(1)=dnull
        ys(1)=dnull
        xs(2)=x
        ys(2)=y
        xmin=dnull
        xmax=dnull
        ymin=dnull
        ymax=dnull
      endif
      do 1 j=1,ngg
      db=w*(yk*x-xk*y-d)
      b=b+db+fekru
      bq=db*db
      e=e+bq
      x=x+dcos(b)*gl
      y=y+dsin(b)*gl
      rmin=dmin1(rmin,bq)
      rmax=dmax1(rmax,bq)
      if(fin) then
        xs(j+2)=x
        ys(j+2)=y
        xmin=min(xmin,x)
        xmax=max(xmax,x)
        ymin=min(ymin,y)
        ymax=max(ymax,y)
      endif
    1 continue
      db=w*(yk*x-xk*y-d)
      b=b+db+fekru
      bq=db*db
      e=e+bq
      xx=x+halb*dcos(b)*gl
      yy=y+halb*dsin(b)*gl
      bb = b/w*gl/10.d0
      ee = e/gl
        rmin=dmin1(rmin,bq)
        rmax=dmax1(rmax,bq)
      if(fin) then
        xs(ngg+3)=xx
        ys(ngg+3)=yy
        xmin=min(xmin,xx)
        xmax=max(xmax,xx)
        ymin=min(ymin,yy)
        ymax=max(ymax,yy)
      endif
      ramin=gl/dsqrt(rmin)
      ramax=gl/dsqrt(rmax)
      return
      end