SUBROUTINE ABQMAIN
C
C	PROGRAM ANALYTICAL
C============================================================================
	PARAMETER (NN=150000)
	DOUBLE PRECISION t(NN),tau(NN),csi(NN),u(NN),theta(NN),gamdot(NN),csai(NN),csaidot(NN)
	DOUBLE PRECISION gamma(NN),x(NN+1),y(NN+1),ksquare,k1square,Omega1,Omega3
        DOUBLE PRECISION coef,u1,u2,u3,beta,h,b,a,I11,I22,I33
C================================================================================
C= RIGID BODY MOTION 3.2.16; Subsection FORCED MOTION OF A RIGID BODY; Case1   ==
C= The roots of f(u),and the elliptic integrals (Omega1e and Omega2e)  are     ==
C= assumed computed as function of I11,I22,I33,omega1,omega2,omega3,theta0,Mg,l==
C= I11=5=I22,I33=1,omega1=0.5,omega2=0.0,omega3=50,theta0=20.0,Mg=20,l=1.0     ==
C= For the case 2, omega1=0.0,omega2=0.0,omega3=50 and the data will change    ==
C================================================================================
	pi=3.141592653589
	Omega1e=1.5722828033360230
	Omega3e=4.168
	beta=8.0
	u3=0.91122001192134283
	u2=0.95172719682607865
	u1=11.607995412038493
	b=9.3969262078590852
	a=10.0
	I11=5.0
	I22=5.0
	I33=1.0
	ksquare=(u2-u3)/(u1-u3)
	k1square=(u1-u2)/(u1-u3)
	Omega1=Omega1e/sqrt(u1-u3)
	Omega3=Omega3e/sqrt(u1-u3)
	h=exp(-Omega3/Omega1*pi)
        coef=4.0*sqrt(h)*sqrt((u1-u3)*(u2-u3))
	csai(1)=0.0*pi/180.0
	gamma(1)=-90.0*pi/180.0        
	dt=0.00010
	DO 20 I=1,NN,1
	 t(I)=(float(I)-1)/10000.0
	 tau(I)=sqrt(beta)/2.0*t(I)-Omega3
         csi(I)=pi*tau(I)/(2.0*Omega1)
	 u(I)=u3+coef*(sin(csi(I))-(h**2)*sin(3.0*csi(I))+(h**6)*sin(5.0*csi(I))-(h**12)*sin(7.0*csi(I)))**2
	 u(I)=u(I)/(1.0-2.0*h*cos(2.0*csi(I))+2.0*(h**4)*cos(4.0*csi(I))-2.0*(h**9)*cos(6.0*csi(I))+2.0*(h**16)*cos(8.0*csi(I)))**2
	 theta(I)=acos(u(I))*180.0/pi
         csaidot(I)=I11/I33*a-u(I)*(b-a*u(I))/(1.0-u(I)**2)
	 csai(I+1)=csai(I)+dt*csaidot(I)
	 gamdot(I)=(b-a*u(I))/(1.0-u(I)**2)
	 gamma(I+1)=gamma(I)+dt*gamdot(I)
	 x(I)=sin(theta(I)*pi/180.0)*cos(gamma(I))
	 y(I)=sin(theta(I)*pi/180.0)*sin(gamma(I)) 
 20	CONTINUE
C============================================================================
C Write the output
C============================================================================
	OPEN(UNIT=98,FILE='angle.out')
	OPEN(UNIT=99,FILE='locus.out')
       	DO 2000 I=1,NN,200
	WRITE(98,1998)t(I),theta(I)
       	WRITE(99,1999)x(I),y(I)
1998	FORMAT(E15.10,' ',E15.10)
1999	FORMAT(E25.10,' ',E25.10)
2000	CONTINUE
	CLOSE(UNIT=98)
	CLOSE(UNIT=99)
C============================================================================
	STOP
	END
C ************* EOF