   	program orbit6n
c       modifizierte Version fuer Perl Programm 8.2.2001
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c       This program calculates the orbit of a test body 
c       (for instance a globular cluster) in the galactic potential
c	of Allen & Santillan (see Rev.Mex.Astr.Astrofis. 22, 255)
c       which models the gravitational potential of the Milky Way.
c
c	Version 6: This version integrates the complete set of 
c       cylindrical coordinates r, z, phi and puts out positions and 
c       velocities in equidistant (!!) time steps.
c
c       Coordinate systems: 
c       External: right-handed galactocentric cartesian coordinates 
c       x,y,z (x from GC to antisolar direction, z to Galactic North Pole) 
c       and associated velocity component U,V,W.
c       Internal: Corresponding cylindricals coordinates r,z,phi
c              
c	Units: kpc(length), km/s(velocity), km**2/s**2(energy);
c	km/(s**2*kpc) (acceleration), 9.778E8y (time, only internal),
c       Gyr (time, external).
c
c       Input: The program requires an input file, which contains 
c       one line with six values each, x,y,z,U,V,W  
c
c       Output: 
c
c       1)
c       Output von t,x,y,z,|R|,rho,U,V,W,v_rho,Theta
c       in das file dummy.txt, was vom Perl-Programm obs2galmulti.pl
c       weiterbearbeitet wird
c
c       2) R_min,R_max,|z|_max,ecc,nze,E_i,jz_i,Theta_i,x_i,y_i
c          z_i,U_i,V_i,W_i nach Standardout
c          wird von Perl in Vergleichsfile von allen Sternen umgeleitet
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c	MAIN ROUTINE
c
	implicit real*8 (a-h,o-z)
	real*8 jz,nze,jz_i,jz_f,th_i
	dimension Y(5)
	character*20 inpfile
	external DERIV
	common/block/jz
c
c       Implizit als real*8 definiert sind:
c       xa,ya,za,vxa,vya,vza,phia,vphia,Tfin,Ti,dt
c       RRmin,RRmax,ecc,zmax,rhozm
c       xx,yy,rho,vx,vy,S,vrho,H,RR,En,TH 

c----------------------------------------------------
c        write(*,*)'       program orbit6n          '
c        write(*,*)' integration of galactic orbits'
c----------------------------------------------------
c       Einlesen der von obs2galn.f generierten Daten 
        read(*,*) xa,ya,za,vxa,vya,vza 
c--------------------------------------------------------
c	Einlesen von Zeit und Zeitschritt aus obs2galmulti.pl
	read(*,*) Tfin
	Tfin=Tfin/0.9778d0

	read(*,*) dt
	dt=dabs(dt)
	dt=dt/0.9778d0
c-----------------------------------------------------------
c
	pi=4.d0*datan(1.d0)
	rho=dsqrt(xa**2+ya**2)
	vrho=(xa*vxa+ya*vya)/rho
	phia=datan2(ya,xa)
	Y(1)=rho
	Y(2)=vrho
	Y(3)=za
	Y(4)=vza
	Y(5)=phia
c--------------------------------------------------------
c       Initialisierung von R_min,max,|z|_max
	RRmin=dsqrt(Y(1)**2+Y(3)**2)
	RRmax=dsqrt(Y(1)**2+Y(3)**2)
        zmax=dabs(za)
        rhozm=rho

c-------z-component of angular momentum vector:
c
	jz=xa*vya-ya*vxa
c
c-------Energy check:
c       Energie und J_z zu Beginn
	call E(rho,vrho,za,vza,En)

        E_i=En
        jz_i=jz
        th_i=-jz/Y(1)
         
c	write(*,*)
c	write(*,*)'angular momentum at initial point:', jz_i
c	write(*,*)'energy at initial point:', En_i
c---------------------------------------------------
 	t=0.d0
	ti=0.d0
	if(Tfin.lt.0.d0) dt=-dt
	H=dt
	m=0
c-------------------------------------------------
c-------Orbit integration loop
c
 100	if(dabs(t).ge.dabs(ti)) goto 300
 120	if(dabs(H).ge.dabs(t-ti)) H=ti-t
c
 	if(dabs(t).ge.dabs(Tfin)) goto 350
	if(dabs(H).ge.dabs(t-Tfin)) H=Tfin-t
c
	call DIFSYA(5,DERIV,1.d-11,H,t,Y)
	if(H.eq.0.d0) then
	write(*,*)'cut off because of H=0.'
	goto 400
	end if
	goto 100
c
c-------------------------------------------------------------------
c-------------------------------------------------------------------
c-------Einige Groessen berechnet
 300 	xx=Y(1)*dcos(Y(5))
	yy=Y(1)*dsin(Y(5))
	TH=-jz/Y(1)
	RR=dsqrt(Y(1)**2+Y(3)**2)
        vx=(Y(2)*xx-jz*yy/Y(1))/Y(1)
	vy=(Y(2)*yy+jz*xx/Y(1))/Y(1)
c-------------------------------------------------------------------
c       Berechnung von Rmin,Rmax,|z|_max,rhozm
        if(RRmin.ge.RR) then
        RRmin=RR
        end if

        if(RRmax.le.RR) then
        RRmax=RR
        endif

        if(zmax.le.(dabs(Y(3)))) then 
        zmax=dabs(Y(3))
        rhozm=Y(1)
        endif
c------------------------------------------------------------------        

c       Output von t,x,y,z,|R|,rho,U,V,W,v_rho,Theta
c       in das file dummy.txt, was vom Perl-Programm obs2galmulti.pl
c       weiterbearbeitet wird
c        open (unit=4,file='dummy.txt',status='unknown') 
c        write(4,310) t*0.9778d0,xx,yy,Y(3),RR,Y(1),vx,vy,Y(4),Y(2),TH
        write(*,310) t*0.9778d0,xx,yy,Y(3),RR,Y(1),vx,vy,Y(4),Y(2),TH

 310	format('b',f10.4,1x,5(f10.6,1x),5(1x,f10.4))

	m=m+1
	ti=m*dt
	goto 120

c        close (unit=4)
c------------------------------------------------------------------
 350	xx=Y(1)*dcos(Y(5))
	yy=Y(1)*dsin(Y(5))
        vx=(Y(2)*xx-jz*yy/Y(1))/Y(1)
	vy=(Y(2)*yy+jz*xx/Y(1))/Y(1)
	S=(Y(5)-phia)/(2.d0*pi)

c-------Final energy check
c
	call E(Y(1),Y(2),Y(3),Y(4),En)

c	write(*,360)t*0.9778d0
c 360	format(' final point of orbit reached at t =',f10.3 )
c	write(*,370) S
c 370	format(' number of revolutions:',f10.2)
        
        E_f=En
        jz_f= xa*vya-ya*vxa

c	write(*,*)'energy at final point:', E_f
c	write(*,*)'angular momentum at final point',jz_f
c------------------------------------------------------------------
c       Ausgabe von R_min,R_max,|z|_max,ecc,nze,E_i,jz_i,S
c       S ist die Zahl der Umlaeufe
        ecc=(RRmax-RRmin)/(RRmax+RRmin)
        nze=zmax/rhozm
          
        write(*,401) RRmin,RRmax,zmax,ecc,nze,E_i,jz_i,S,th_i,
     $  xa,ya,za,vxa,vya,vza
 401	format('a',5(f10.4,1x),1(e10.4,1x),3(f10.4,1x),1x,6(f10.4,1x))
c---------------------------------------------------------------------
c-------------------------------------------------------------------
c
 400	stop

	end
c-----------------------------------------------------------------
c----------------------------------------------------------------
c---------------------------------------------------------------
c--------------------------------------------------------------
c-------------------------------------------------------------------
c
	SUBROUTINE DERIV(t,Y,DYDt)
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c       This subroutine builds up the equations of motion
c       (for the forces see Allen&Santillan, Rev.Mex. 22, 255-263 )
c
c	Implementation 14.2.1995 , M.Odenkirchen
c
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
	implicit real*8 (a-h,o-z)
	real*8 jz,M1,M2,M3,MinR
	common/block/jz
	dimension Y(5),DYDt(5)
c
	r=Y(1)
	z=Y(3)
c
c	Partial Derivatives of the Potential:
c
c	Component 1 (central mass):
	M1=606.0d0
	b1=0.3873d0
	root1=dsqrt(r*r+z*z+b1*b1)
	root13=root1*root1*root1
	dp1dr=r*M1/root13
	dp1dz=z*M1/root13
c
c	Component 2 ("disk"):
	M2=3690.d0
	a2=5.3178d0
	b2=0.2500d0
	root2=dsqrt(z*z+b2*b2)
	root3=dsqrt(r*r+(a2+root2)**2)
	root3=root3**3
	dp2dr=r*M2/root3
	dp2dz=z*(M2/root3)*((a2+root2)/root2)
c
c	Component 3 (halo):
	M3=4615.0d0
	a3=12.0d0
	rr=dsqrt(r*r+z*z)
	rr3=rr*rr*rr
	x=rr/a3
	MinR=M3*(x**2.02)/(1d0+x**1.02)
	dp3dr=r*MinR/rr3
	dp3dz=z*MinR/rr3
c
c
	DYDt(1)=Y(2)
	DYDt(2)=100d0*(-dp1dr-dp2dr-dp3dr)+jz**2/r**3
	DYDt(3)=Y(4)
	DYDt(4)=(-dp1dz-dp2dz-dp3dz)*100d0
	DYDt(5)=jz/Y(1)**2
c
	return
	end
c
c
c
c
c
c
c
	SUBROUTINE E(r,vr,z,vz,En)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c	Specific total energy En of a particle at (r,z) 
c       with velocity (vr,vz) and z-angular momentum jz 
c	in the potential of Allen&Santillan (Rev.Mex. 22, 255)
c
c	Implementation 13.2.1995, M.Odenkirchen
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
	implicit real*8 (a-h,o-z)
	real*8 jz,M1,M2,M3,MinR
	common/block/jz
c
	r2=r*r
	z2=z*z
	rr=dsqrt(r2+z2)
c
c	Component 1 (central mass):
	M1=606.0d0
	b1=0.3873d0
	root1=dsqrt(r2+z2+b1*b1)
	p1=-M1/root1
c
c	Component 2 ("disk"):
	M2=3690.d0
	a2=5.3178d0
	b2=0.2500d0
	root2=dsqrt(z2+b2*b2)
	root3=dsqrt(r2+(a2+root2)**2)
	p2=-M2/root3
c
c	Component 3 (halo):
	M3=4615.0d0
	a3=12.0
	rr=dsqrt(r2+z2)
	x=rr/a3
	xx=1d0+x**1.02
	MinR=M3*x**2.02/xx
	p3=-MinR/rr
c
	const=m3/(1.02*a3)
	brac1=-(1.02/xx)+dlog(xx)
	xx2=1d0+(100d0/a3)**1.02
	brac2=-(1.02/xx2)+dlog(xx2)
	p3=p3-const*(brac2-brac1)
c
c	Total potential:
	p=100d0*(p1+p2+p3)
c
c	factor 100 for adjustment of the potential
c	to unit km**2/s**2 !!
c
c	Total energy:
	En=.5d0*(vr**2+(jz/r)**2+vz**2)+p
c
	return
	end
c
c
c
      SUBROUTINE DIFSYA (N,F,EPS,H,X,Y)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c       A modified version of the Bulirsch-Stoer integrator
c       to be used for orbit integration
c
c       M.Odenkirchen
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      real*8 Y(N),YA(06),YL(06),YM(06),DY(06),DZ(06),DT(06,7),
     1D(7),S(06),X,XN,H,G,B,B1,U,V,C,TA,W,DABS,DFLOAT,EPS
      dimension EP(4)
      logical   KONV,BO,KL,GR
      data EP/0.4E-1,0.16E-2,0.64E-4,0.256E-5/
c
      JTI=0
      FY=1.
      DO 100 I=1,N
  100 YA(I)=Y(I)
      CALL F (X,Y,DZ)
   10 XN=X+H
      BO=.FALSE.
      DO 110 I=1,N
  110 S(I)=0.D0
      M=1
      JR=2
      JS=3
      DO 260 J=1,10
      IF(.NOT.BO) GOTO 200
      D(2)=1.777777777777778D0
      D(4)=7.111111111111111D0
      D(6)=2.844444444444444D1
      GOTO 201
  200 D(2)=2.25D0
      D(4)=9.D0
      D(6)=3.6D1
  201 IF(J.LE.7) GOTO 202
      L=7
      D(7)=6.4D1
      GOTO 203
  202 L=J
      D(L)=M*M
  203 KONV=L.GT.3
      M=M+M
      G=H/DFLOAT(M)
      B=G+G
      M=M-1
      DO 210 I=1,N
      YL(I)=YA(I)
  210 YM(I)=YA(I)+G*DZ(I)
	if (YM(1).lt.0.d0) then
	H=H*0.5d0
	goto 10
	end if
      DO 220 K=1,M
      CALL F (X+DFLOAT(K)*G,YM,DY)
      DO 220 I=1,N
      U=YL(I)+B*DY(I)
      YL(I)=YM(I)
      YM(I)=U
      U=DABS(U)
      IF(U.GT.S(I)) S(I)=U
  220 CONTINUE
	if (YM(1).lt.0.d0) then
	H=H*0.5d0
	goto 10
	end if
      CALL F (XN,YM,DY)
      KL=L.LT.2
      GR=L.GT.5
      FS=0.
      DO 233 I=1,N
      V=DT(I,1)
      C=(YM(I)+YL(I)+G*DY(I))*0.5D0
      DT(I,1)=C
      TA=C
      IF(KL) GOTO 233
      DO 231 K=2,L
      B1=D(K)*V
      B=B1-C
      W=C-V
      U=V
      IF(B.EQ.0.D0) GOTO 230
      B=W/B
      U=C*B
      C=B1*B
  230 V=DT(I,K)
      DT(I,K)=U
  231 TA=U+TA
      IF(.NOT.KONV) GOTO 232
      IF(DABS(Y(I)-TA).GT.EPS) KONV=.FALSE.
  232 IF(GR.OR.S(I).EQ.0.D0) GOTO 233
      FV=DABS(W)/S(I)
      IF(FS.LT.FV) FS=FV
  233 Y(I)=TA
      IF(FS.EQ.0.D0) GOTO 250
      FA=FY
      K=L-1
      FY=(EP(K)/FS)**(1./FLOAT(L+K))
      IF(L.EQ.2) GOTO 240
      IF(FY.LT.0.7*FA) GOTO 250
  240 IF(FY.GT.0.7) GOTO 250
      H=H*FY
      JTI=JTI+1
      IF(JTI.GT.5) GOTO 30
      GOTO 10
  250 IF(KONV) GOTO 20
      D(3)=4.D0
      D(5)=1.6D1
      BO=.NOT.BO
      M=JR
      JR=JS
  260 JS=M+M
      H=H*0.5D0
      GOTO 10
   20 X=XN
      H=H*FY
      RETURN
   30 H=0.D0
      DO 300 I=1,N
  300 Y(I)=YA(I)
      RETURN
        end


