c     geometry of astatic leaf-spring suspensions   nov. 2000 ew
c
c     Comments and variable names translated by
c     Brett Nordgren – January 3, 2008
c
      program leaffind
c
c     internal parameters:
c     acc: required accuracy for the iteration (cm)
c     xinc:  increment when calculating the partial derivatives (cm)
c     hstep: angle increment for calculating the potential energy (rad)
c     maxit:  abort the iteration after 'maxit' calls to blf
c
c     note: Three different normalizations are used for the rotational moments.
c     Routine ‘blf’ computes with dimensionless torque. With the
c     inversion under the starting parameters this is divided by 5.
c     The physical torque in kp*cm arises as a result of multi-
c     plication by the factor 2*pi/180. The same factor is used with the
c     applied force.
c     Similarly the force also appears with two different normalizations.
c
      implicit double precision (a-h,o-z)
      logical circle
      dimension gg(3),hh(6),f(2,3,2),a(3),b(3),c(3)
      common acc,xinc,ax,pi4,half,ngg,gl,ramin,ramax,fl,itera,hstep,
      -nix,bu,maxit,fakt,dfak,fekru,torad,dnull,itotal,itrace

      itrace=0
      strmax=280.d0
c  1/strmax is the maximum permitted strain of the spring material
      dnull=0.d0
      half=5.d-1
      pi4=datan2(half,half)
c   pi4 = datan2(x,y) = arg(x+yi) = arg(0.5+0.5i) = pi/4
      torad=pi4/45.d0
c   torad = pi/180
      acc=1.d-10
      xinc=1.d-4
      maxit=333
      nix=0
      hstep=1.d-2
      step=1.d-2
      ng=64
      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 ==> leaffind ==> leaf.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 circular spring (xo,yo ignored)')
c    &    /'      = 4, compute partial derivatives of f^2 and mass'/)

   10 read (7,1,err=50,end=50) iopt,xu,yu,wu,xo,yo,wo,fl,ax,fekru,
     & width,thick,emodul,xki,yki,di
    1 format(i1,1x,15f5.3)
      if(iopt.eq.0) then
      write(6,150)
  150 format(' Option iopt=0 or blank in col.1 - stop')
      stop
      endif
      if(iopt.le.4) goto 70
   50 if(iopt.eq.0) then
      write(6,150)
      stop
      endif
      write(6,60)
   60 format(' Input error. The file leaf.in must contain the following'
     &/' 16 parameters in one line in format (i1,1x,15f5.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,width,thick,emodul,
     & xki,yki,di
      write(8,4)     xu,yu,wu,xo,yo,wo,fl,ax,fekru,width,thick,emodul,
     & xki,yki,di
    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/)

      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'/)

      raref=thick*strmax/2.d0
      fakt=emodul*width*thick**3/24.d0
c   fakt = EI/2
      dfak=fakt*2.d0*torad
c   dfak = EI * pi/180 
      gl=fl/ng

      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 ignored
     &d'/' and replaced by values calculated from the spring geometry.'
     &/' The same applies to xki, yki, and di.'/)
      endif

      write(6,2) width,thick*10.d0,emodul/100.d0
      write(8,2) width,thick*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.or.iopt.eq.4) goto 110
      write(6,*) 'Iterative search:'
      write(8,*) 'Iterative search:'
  302 nix=0

   90 format(/' itera      fak        abold        abnew         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
      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)

      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 at
     & 5 cm from hinge.'/)
      if(iopt.eq.3) then
      xo=xu+sehn
      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)

      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)')
      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 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,half,ngg,gl,ramin,ramax,fl,itera,hstep,
      -nix,bu,maxit,fakt,dfak,fekru,torad,dnull,itotal,itrace
c   note... angle unit for b is 10 degrees
c   note... true rotational moment is 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
      zpiq = zpiq*zpiq
      rm=5.d0
      ter=dnull
c   compute the elastic potential in 5 neighbouring positions
      do 30 l=1,5
      h = hstep*(3-l)

      call circle(xo,yo,b,h,xz,yz,bz,r)

      call search(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   rotational moment and seismic mass
      dreh=(e(4)-e(2))/(2.*hstep)
      rmg=dreh/dcos(ax)
      gm=rmg/rm
      hh(5)=gm
c   add the gravity potential
 
      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   rotational moment as derivative of the potentials w.r.t. the rotation angle
      do 70 l=1,4
   70 fkr(l) =(e(l+1)-e(l))/hstep
c   restoring moment as derivative of the rotational moments
      do 80 l=1,3
      fko(l) =(fkr(l+1)-fkr(l))/hstep
c   square of the natural frequency
   80 fq(l) = fko(l)/gm*981./rm**2/zpiq
      fm=fq(2)
      fdif=0.5d0*(fq(3)-fq(1))/(hstep*rm)
      ter=(fq(3)-fq(2)*2.+fq(1))/(hstep*rm)**2
      return

   40 fm=dnull
      fdif=dnull
      50 j=1,3
      gg(j)=dnull
      hh(j+3)=dnull
   50 hh(j)=dnull
      return
      end
      subroutine circle(xo,yo,ba,defl,x,y,b,r)
      implicit double precision (a-h,o-z)
      common acc,xinc,ax,pi4,half,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 search(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,half,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
      abnew=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 xkold=xk
      ykold=yk
      dold=d
      abold=abnew
      xold=x
      yold=y
      bold=b
    3 call matin(a,dif,abnew*fak,xkor)
      xk=xkold+xkor(1)
      yk=ykold+xkor(2)
      d= dold+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
      abnew=ab(dif(1),dif(2),dif(3))
      if(itrace.eq.1) write(6,4) itera,fak,abold,abnew,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(abnew.le.abq) return
      if(abnew.lt.abold) goto 2
      if(itera.gt.2) goto 14
      nix=0
      goto 7

      14 fak=fak*32.d00
      dif(1)=xz-xold
      dif(2)=yz-yold
      dif(3)=bz-bold
      goto 3
      end
      subroutine matin(a,dif,wq,xkor)
    c   matrix invert?
      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)
    c   Gaussian elimination?
      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)
    c   Balanced loss function???
      implicit double precision (a-h,o-z)
      common acc,xinc,ax,pi4,half,ngg,gl,ramin,ramax,fl,itera,hstep,
      -nix,bu,maxit,fakt,dfak,fekru,torad,dnull,itotal,itrace
      itera=itera+1
      xk = xkk
      yk = ykk
      d = dk*5.d0
      w=torad*gl
      e=dnull
      b=10.d0*bu*w/gl
      x=half*dcos(b)*gl
      y=half*dsin(b)*gl
      rmin=1.d12
      rmax=dnull
      do 1 j=1,ngg
      db=w*(yk*x-xk*y-d)
      b=b+db+fekru
      bq=db*db
      e=e+bq
      rmin=dmin1(rmin,bq)
      rmax=dmax1(rmax,bq)
      x=x+dcos(b)*gl
    1 y=y+dsin(b)*gl
      db=w*(yk*x-xk*y-d)
      b=b+db+fekru
      bq=db*db
      e=e+bq
      rmin=dmin1(rmin,bq)
      rmax=dmax1(rmax,bq)
      xx=x+half*dcos(b)*gl
      yy=y+half*dsin(b)*gl
      bb = b/w*gl/10.d0
      ee = e/gl
      ramin=gl/dsqrt(rmin)
      ramax=gl/dsqrt(rmax)
      return
      end