Fortran Code for The Trebuchet

c	MonteCarloTreb1.f
c	Calculates the range of a trebuchet with a free sling
c	and a hinged counterweight.  Assumes random input parameters
c	for a monte-carlo exploration of the possible designs.
c	Does stage 1 of the manuscript, 5 or so throws per second.
c
c	D. B. Siano             dimona@home.net       July 9, 1998
c	Documentation is at http://www.members.home.com/dimona
c	
c	Written for Absoft Fortran 77
c 	Calculates the second derivatives of the angles using
c	expressions derived from a Mathematica program.
c	Uses fourth-order Runge Kutta for second order DE
c	formula 25.5.20 in Abromowitz's "Handbook of Mathemtical
c	Functions" NBS, reprinted by Dover (1968).
c	Utilizes a variable stepsize (set by an accuracy tolerance, eps)
	Real*8 l1,l2,l3,l4,l5,m1,m2,mb,g,r,vx,vy
	Real*8 frth,tmax,thmax,phimax,psimax,rmax,rmaxth
	Real*8 ths,phis,psis,thds,phids,psids,hz,h1,eps
	Real*8 psiz,thz,phiz,thdz,phidz,psidz,delth
	Real*8 hstep,phist2,phidst2,psist2,psidst2
	Real*8 thst,thst2,thdst2,phist,phidst,psist,psidst
	Real*8 der1,der2,der3,thdst,h,psizdeg
	Real*8 a,b,rat,pecw,keproj,eneff,vxmax,vymax,thgo,effmax
	COMMON l1,l2,l3,l4,m1,m2,mb,g
		
c	Input the data here--uses ft, pounds
	l1=1.33
	l2=5.59
	l3=5.59
	l4=2.49
	l5=4.84
	m1=100.
	m2=1.
	mb=5.

c	write out the input data
	write(*,904)
904	Format(' free sling, hinged cw model')
	write(*,907)
907	Format(' MonteCarloTreb1.f')

c	acceleration of gravity
	g=32.
c	tolerance factor:  increase this for greater speed, less accuracy
	eps=1.d-7
	
	write(*,908)
908	Format(' j,l1,l2,l3,l4,l5,m1,m2,mb,rmax,eff,eneff,thgo,
     - vymax,vxmax,thmaxdeg,phimaxdeg,psimaxdeg,ncount,tmax,h,psizdeg')
	effmax=0.
	
	psiz=dasin(l5/l2)
c	set mone=1 to do a single throw, using the input values above
c	set mone=2 to do multiple random throws, using the input values
c	calculated below
	mone=2
	if(mone.eq.2) go to 51
	If(mone.eq.1) go to 50
	
51	continue
c	set the upper and lower limits on the values of psiz that are to be
c	used for the random variation here
	a=3.1415926/6
	b=3.1415926/3
	
c	set the number of throws here	
	Do 101 j=1,200


c	choose psiz between 30 and 60
	k=5*j
	psiz=(b-a)*RAN2(k)+a
	
c	assume cw falls two ft
	l1=2./(1.+dsin(psiz))
	
c	assume l2/l1 is between 3 and 5
	l=j
	rat=(5.-3.)*RAN2(l)+3.
	l2=rat*l1
	l5=l2*dsin(psiz)


c	assume sling length is between .5*l2 and 1.5*l2
	m=7*j
	l3=((1.5-.5)*RAN2(m)+.5 )*l2
	
c	choose l4 to be between .5*l1 and 1.*l1
	mm=9*j
	l4=((1.-.5)*RAN2(mm)+.5)*l1
	
50	continue


c	the starting angles: cw hangs, sling is horizontal
	thz=psiz+3.1415926/2.
	phiz=(3.1415926/2)-psiz
	
c	the angles start with zero rate of change
c	the "d" means first derivative, "dd" means the second derivative
	thdz=0.
	phidz=0.
	psidz=0.
	
c	Calculate pts from t=0 to t=hz seconds	
c	with a starting assumed step size of hz/hstep
c	you may need to adjust hz for weird cases
	hz=2.*dsqrt(2.*l2/g)
	hstep=600.
	h1=hz/hstep
	
c	write(*,988) hz, eps
988	format(' assumed max time for throw =',F8.2,', eps=',E8.1)
c	write(*,909) thz,phiz,psiz
909	format(' strt angles=',3d14.4)

c	der1,2,3 are used for checking correctness of dd functions
	der1=thdd(thz,thdz,phiz,phidz,psiz,psidz)
	der2=phidd(thz,thdz,phiz,phidz,psiz,psidz)
	der3=psidd(thz,thdz,phiz,phidz,psiz,psidz)
c	write(*,906) der1,der2,der3
906	format('der1,2,3=',3d14.4)

	rmax=0
c	the theoretical maximum range
	rmaxth=2*m1*(1.+dsin(psiz))*l1/m2
	tim=0.
c	save start vals
	ths=thz
	phis=phiz
	psis=psiz
	thds=thdz
	phids=phidz
	psids=psidz
	
c	counts the number of steps taken
	ncount=0
	
c	Primary loop:
c	get the maximum range for this set of parameters
c	use ncount to avoid infinite loops
c	the runge-kutta variable step size algorithm
	do while(tim < hz.AND.ncount < 4000)
		ncount=ncount+1
		tim=tim+2*h1
c		1 step of size 2*h1:	
	Call thstep(2*h1,ths,thds,phis,phids,psis,psids,thst2,thdst2)
	Call phistep(2*h1,ths,thds,phis,phids,psis,psids,phist2,phidst2)
	Call psistep(2*h1,ths,thds,phis,phids,psis,psids,psist2,psidst2)

c	now do 2 steps of size h1:
	do 100 k=1,2
	  Call thstep(h1,ths,thds,phis,phids,psis,psids,thst,thdst)	
	  Call phistep(h1,ths,thds,phis,phids,psis,psids,phist,phidst)
	  Call psistep(h1,ths,thds,phis,phids,psis,psids,psist,psidst)
	  
	  ths=thst
	  thds=thdst
	  phis=phist
	  phids=phidst
	  psis=psist
	  psids=psidst
	  delth=thst2-thst
	  frth=delth/thst
100	continue

c	Only consider the first rotation of the beam till it is
c	45 deg past vertical
	If(thst < -.707) goto 888
c	Only consider solutions with psimaxdeg<270:no twirling
	If(psist > 4.71) goto 888
	
c	Compare the two results and adjust the stepsize accordingly	
	if(dabs(frth).LT.eps)  then
	h1=h1*2
	else 
	h1=h1/2
	endif
	
c	calculate the range assuming release now
	vx=-(l3*(psidst-thdst)*dcos(psist-thst))-l2*thdst*dcos(thst)
    	vy= l3*(psidst-thdst)*dsin(psist-thst)-l2*thdst*dsin(thst)
     	r=2*vx*vy/g
	
     	if(r > rmax) then
	
68	format(' range=',2d14.4)
		rmax=r
		tmax=tim
		thmax=thst2
		phimax=phist2
		psimax=psist2
		vxmax=vx
		vymax=vy
	endif

	end do
888	continue
c	calculate the energetics, other characteristics of the throw
	keproj=m2*(vxmax*vxmax+vymax*vymax)/2.
	pecw=m1*g*l1*(1.+Sin(psiz))
	eneff=keproj/pecw
	thmaxdeg=thmax*57.29578
	phimaxdeg=phimax*57.29578
	psimaxdeg=psimax*57.29578
	eff=rmax/rmaxth
	vyy=vymax
	yxx=vxmax
	thgo=57.29578*datan(vymax/vxmax)
c	make the distance of cw fall, h, a constant
	h=l1*(1.+ sin(psiz))
	
c	look only for efficient trebs
	if(eff > 0) then
	effmax=eff
	psizdeg=psiz*57.29578
	write(*,900) j,l1,l2,l3,l4,l5,m1,m2,mb,rmax,eff,eneff,thgo,
     - vymax,vxmax,thmaxdeg,phimaxdeg,psimaxdeg,ncount,tmax,h,psizdeg
	endif
	
900	Format(I4,5(F6.2),F6.1,2(F5.1),f7.1,2f6.3,3(f6.1),3(f6.0),
     -	I5,f6.3,2(f6.1))
903	Format(' start angles =',3F8.3,' rad')
902	Format(//,'input paras l,m=',5F7.2,3F7.1)		
c	8	continue
c	9	Continue
101	continue
	pause
	stop
	end
	
c-----------------------------------------------------------------
c	random number generator	
         FUNCTION RAN2(IDUM)
         PARAMETER (M=714025,IA=1366,IC=150889,RM=1.4005112E-6)
         DIMENSION IR(97)
         DATA IFF /0/
         IF(IDUM.LT.0.OR.IFF.EQ.0)THEN
         IFF=1
         IDUM=MOD(IC-IDUM,M)
         DO 11 J=1,97
           IDUM=MOD(IA*IDUM+IC,M)
           IR(J)=IDUM
11       CONTINUE
         IDUM=MOD(IA*IDUM+IC,M)
         IY=IDUM
         ENDIF
         J=1+(97*IY)/M
         IF(J.GT.97.OR.J.LT.1)PAUSE
         IY=IR(J)
         RAN2=IY*RM
         IDUM=MOD(IA*IDUM+IC,M)
         IR(J)=IDUM
         RETURN
         END
c-------------------------------------------------------------------
c	the three second derivatives of the angles as functions	     
	function thdd(th,thd,phi,phid,psi,psid)
	Real*8 l1,l2,l3,l4,m1,m2,mb,th,thd,phi,phid,psi,psid,g
	Real*8 thdd
	common l1,l2,l3,l4,m1,m2,mb,g
	dsphi=dsin(phi)
	dcphi=dcos(phi)
	dspsi=dsin(psi)
	dcpsi=dcos(psi)
	dsth=dsin(th)
	dcth=dcos(th)
            thdd=(-3*(-2*l1*l4*m1*phid**2*dsphi - 
     -      4*l1*l4*m1*phid*thd*dsphi - 
     -      2*l1*l4*m1*thd**2*dsphi + 
     -      2*l1**2*m1*thd**2*dcphi*dsphi + 
     -      2*l2*l3*m2*psid**2*dspsi - 
     -      4*l2*l3*m2*psid*thd*dspsi + 
     -      2*l2*l3*m2*thd**2*dspsi - 
     -      2*l2**2*m2*thd**2*dcpsi*dspsi + 
     -      2*g*l2*m2*dcpsi*(dspsi*dcth-dcpsi*dsth) - 
     -      2*g*l1*m1*dsth + 2*g*l2*m2*dsth - 
     -      g*l1*mb*dsth + g*l2*mb*dsth + 
     -      2*g*l1*m1*dcphi*(dsphi*dcth+
     -      dcphi*dsth)))/
     -  (2.*(-3*l1**2*m1 - 3*l2**2*m2 - l1**2*mb + 
     -      l1*l2*mb - l2**2*mb + 3*l1**2*m1*dcphi**2 + 
     -      3*l2**2*m2*dcpsi**2))
  	 return
  	 end
	 
   	function phidd(th,thd,phi,phid,psi,psid)
   	Real*8 l1,l2,l3,l4,m1,m2,mb,th,thd,phi,phid,psi,psid,g
	Real*8 phidd
   	common l1,l2,l3,l4,m1,m2,mb,g
   	dsphi=dsin(phi)
	dcphi=dcos(phi)
	dspsi=dsin(psi)
	dcpsi=dcos(psi)
	dsth=dsin(th)
	dcth=dcos(th)
         phidd=-((-(l1*thd**2*dsphi) - g*(dsphi*dcth
     -    +dcphi*dsth))/l4) + 
     -  (3*(l4 - l1*dcphi)*
     -     (-2*l1*l4*m1*phid**2*dsphi - 
     -       4*l1*l4*m1*phid*thd*dsphi - 
     -       2*l1*l4*m1*thd**2*dsphi + 
     -       2*l1**2*m1*thd**2*dcphi*dsphi + 
     -       2*l2*l3*m2*psid**2*dspsi - 
     -       4*l2*l3*m2*psid*thd*dspsi + 
     -       2*l2*l3*m2*thd**2*dspsi - 
     -       2*l2**2*m2*thd**2*dcpsi*dspsi + 
     -       2*g*l2*m2*dcpsi*(dspsi*dcth-dcpsi*dsth) - 
     -       2*g*l1*m1*dsth + 2*g*l2*m2*dsth - 
     -       g*l1*mb*dsth + g*l2*mb*dsth + 
     -       2*g*l1*m1*dcphi*(dsphi*dcth+
     -       dcphi*dsth)))/
     -   (2.*l4*(-3*l1**2*m1 - 3*l2**2*m2 - l1**2*mb + 
     -       l1*l2*mb - l2**2*mb + 3*l1**2*m1*dcphi**2 + 
     -       3*l2**2*m2*dcpsi**2))
   	return
   	end
   	
   	function psidd(th,thd,phi,phid,psi,psid)
   	Real*8 l1,l2,l3,l4,m1,m2,mb,th,thd,phi,phid,psi,psid,g
	Real*8 psiddtop,psiddbot,psidd
   	common l1,l2,l3,l4,m1,m2,mb,g
   	dsphi=dsin(phi)
	dcphi=dcos(phi)
	dspsi=dsin(psi)
	dcpsi=dcos(psi)
	dsth=dsin(th)
	dcth=dcos(th)
         psiddtop=
     -   l3*m2*(-(l4*m1*(l4 - l1*dcphi)*(l4**2*m1 - 
     -    l1*l4*m1*dcphi)) + 
     -     l4**2*m1*(l1**2*m1 + l4**2*m1 + l2**2*m2 + l3**2*m2 +
     -   (l1**2*mb)/3.
     -       - (l1*l2*mb)/3.+(l2**2*mb)/3.-2*l1*l4*m1*dcphi- 
     -    2*l2*l3*m2*dcpsi))*
     -   (-(l2*thd**2*dspsi) + g*(dspsi*dcth-dcpsi*dsth)) - 
     -  l3*m2*(-l3 + l2*dcpsi)*(-(l4*m1*(l4**2*m1 - 
     -   l1*l4*m1*dcphi)*
     -        (-(l1*thd**2*dsphi) - g*(dsphi*dcth+
     -    dcphi*dsth))) + 
     -     l4**2*m1*(l1*l4*m1*phid**2*dsphi + 
     -    2*l1*l4*m1*phid*thd*dsphi
     -        -l2*l3*m2*psid**2*dspsi +
     -    2*l2*l3*m2*psid*thd*dspsi - 
     -        g*l3*m2*(dspsi*dcth-dcpsi*dsth) + g*l1*m1*dsth 
     -   - g*l2*m2*dsth + 
     -        (g*l1*mb*dsth)/2. - (g*l2*mb*dsth)/2. -
     -      g*l4*m1*(dsth*dcphi+dcth*dsphi)))
     
        psiddbot=-(l1**2*l3**2*l4**2*m1**2*m2) -
     -   l2**2*l3**2*l4**2*m1*m2**2 - 
     -  (l1**2*l3**2*l4**2*m1*m2*mb)/3. + 
     -    (l1*l2*l3**2*l4**2*m1*m2*mb)/3. - 
     -  (l2**2*l3**2*l4**2*m1*m2*mb)/3. + 
     -   l1**2*l3**2*l4**2*m1**2*m2*dcphi**2 + 
     -  l2**2*l3**2*l4**2*m1*m2**2*dcpsi**2
         psidd=psiddtop/psiddbot
         return
         end
	 
c	The Runge-Kutta formulas as functions

        Subroutine thstep(h,th,thd,phi,phid,psi,psid,thst,thdst)
	Real*8 k1,k2,k3,k4,h,th,thd,phi,phid,psi,psid,thst,thdst
        k1=h*thdd(th,thd,phi,phid,psi,psid)
        k2=h*thdd(th+h/2*thd+h*k1/8,thd+k1/2,phi,phid,psi,psid)
        k3=h*thdd(th+h/2*thd+h*k1/8,thd+k2/2,phi,phid,psi,psid)
        k4=h*thdd(th+h*thd+h*k3/2,thd+k3,phi,phid,psi,psid)
	thst=th+h*(thd+(k1+k2+k3)/6)
        thdst=thd+(k1+2*k2+2*k3+k4)/6
        return
        end
      
	Subroutine phistep(h,th,thd,phi,phid,psi,psid,phist,phidst)
	Real*8 k1,k2,k3,k4,h,th,thd,phi,phid,psi,psid,phist,phidst
        k1=h*phidd(th,thd,phi,phid,psi,psid)	
	k2=h*phidd(th,thd,phi+h/2*phid+h*k1/8,phid+k1/2,psi,psid)
	k3=h*phidd(th,thd,phi+h/2*phid+h*k1/8,phid+k2/2,psi,psid)
	k4=h*phidd(th,thd,phi+h*phid+h*k3/2,phid+k3,psi,psid)
	phist=phi+h*(phid+(k1+k2+k3)/6)
	phidst=phid+(k1+2*k2+2*k3+k4)/6
        return
        end

	Subroutine psistep(h,th,thd,phi,phid,psi,psid,psist,psidst)
	Real*8 k1,k2,k3,k4,h,th,thd,phi,phid,psi,psid,psist,psidst
	k1=h*psidd(th,thd,phi,phid,psi,psid)	
	k2=h*psidd(th,thd,phi,phid,psi+h/2*psid+h*k1/8,psid+k1/2)
    	k3=h*psidd(th,thd,phi,phid,psi+h/2*psid+h*k1/8,psid+k2/2)
	k4=h*psidd(th,thd,phi,phid,psi+h*psid+h*k3/2,psid+k3)
	psist=psi+h*(psid+(k1+k2+k3)/6)
	psidst=psid+(k1+2*k2+2*k3+k4)/6
        return
        end

Here is the output from the program:

  free sling, hinged cw model
 MonteCarloTreb1.f
 j,l1,l2,l3,l4,l5,m1,m2,mb,rmax,eff,eneff,thgo, vymax,vxmax,thmaxdeg,phimaxdeg,psimaxdeg,ncount,tmax,h,psizdeg
   1  1.29  5.30  3.73  1.02  2.89 100.0  1.0  5.0  255.4  .639  .640  43.3  62.0  65.9   19.  198.  151.  356  .485   2.0  33.0
   2  1.32  4.99  6.75   .76  2.55 100.0  1.0  5.0  122.5  .306  .341  58.0  56.0  35.0  -21.  217.  111. 1358  .607   2.0  30.7
   3  1.11  3.65  2.94   .87  2.90 100.0  1.0  5.0  210.0  .525  .533  39.9  53.0  63.3    9.  219.  138.  576  .474   2.0  52.6
   4  1.20  5.51  5.47   .61  3.66 100.0  1.0  5.0  238.7  .597  .597  46.2  63.2  60.5   17.  207.  153.  450  .563   2.0  41.6
   5  1.11  3.54  2.98   .82  2.86 100.0  1.0  5.0  210.5  .526  .533  40.6  53.7  62.7    8.  223.  137.  586  .474   2.0  54.1
   6  1.12  4.00  4.58   .67  3.16 100.0  1.0  5.0  254.5  .636  .638  42.9  61.5  66.2    2.  239.  136.  725  .525   2.0  52.3
   7  1.21  4.14  5.95  1.00  2.67 100.0  1.0  5.0  190.8  .477  .478  43.1  53.5  57.1  -28.  256.  103.  752  .558   2.0  40.3
   8  1.24  3.93  3.85   .82  2.43 100.0  1.0  5.0  251.8  .629  .630  43.6  62.0  65.0    5.  227.  137.  601  .469   2.0  38.3
   9  1.18  4.90  5.62   .78  3.42 100.0  1.0  5.0  255.5  .639  .639  45.6  64.6  63.3    4.  223.  138.  603  .556   2.0  44.3
  10  1.22  3.71  2.15  1.09  2.39 100.0  1.0  5.0  149.5  .374  .404  33.9  40.1  59.7   24.  184.  136.  517  .415   2.0  40.1
  11  1.23  4.62  6.65   .81  2.92 100.0  1.0  5.0  179.6  .449  .449  44.8  53.4  53.8  -29.  250.  103. 1058  .596   2.0  39.2
  12  1.20  3.76  3.62  1.01  2.49 100.0  1.0  5.0  267.0  .667  .669  42.8  62.9  67.9    4.  220.  136.  697  .464   2.0  41.4
  13  1.08  4.32  3.03   .76  3.66 100.0  1.0  5.0  165.7  .414  .423  39.3  46.6  56.9   15.  224.  143.  519  .535   2.0  57.9
  14  1.30  6.12  5.58   .92  3.28 100.0  1.0  5.0  244.9  .612  .614  47.3  65.1  60.2   17.  194.  153.  442  .552   2.0  32.4
  15  1.18  5.34  4.75   .80  3.68 100.0  1.0  5.0  264.5  .661  .661  45.1  65.2  64.9   18.  206.  152.  388  .542   2.0  43.5
  16  1.11  3.35  2.43   .83  2.68 100.0  1.0  5.0  162.4  .406  .420  37.5  44.6  58.2   26.  183.  131.  579  .431   2.0  53.2
  17  1.08  3.94  3.82   .80  3.38 100.0  1.0  5.0  246.9  .617  .620  42.4  60.1  65.7    9.  225.  140.  593  .528   2.0  59.0
  18  1.11  3.66  3.36   .74  2.95 100.0  1.0  5.0  234.0  .585  .588  42.2  58.3  64.3    8.  228.  139.  571  .485   2.0  53.6
  19  1.15  5.76  5.92  1.06  4.24 100.0  1.0  5.0  273.7  .684  .685  45.9  67.2  65.2   15.  198.  151.  479  .601   2.0  47.4
  20  1.18  4.46  6.07   .66  3.08 100.0  1.0  5.0  214.6  .537  .537  43.8  57.4  59.9  -16.  249.  117.  921  .572   2.0  43.6

Go Back