	program obs2galn
c       von mir modifizierte Version 30.1.2001
!       Dehnen & Binney + IAU -Werte eingebaut. RN 14/9/05
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c       Conversion of observed kinematic parameters into 
c       galactic parameters to be used in an orbit computation
c
c       Input: (from user-specified file) 
c       equatorial position, distance from Sun, heliocentric 
c       radial velocity and observed absolute proper motions in one line 
c       (format: iah iam as idg idm ds d vr pma pmd
c        with 
c        iah= hours of RA, iam= minutes of RA, as=seconds of RA
c        idg=degrees of DEC (with sign!!), idm= ..., ds=...
c        d= distance in kpc, vr radial velocity in km/s
c        pma= mu_alpha cos delta, pmd = mu_delta  in mas/y )
c	
c       Output: (to file "obs2gal.dat") 
c       galactocentric cartesian coordinates x,y,z (in kpc)
c       and velocity components U,V,W (in km/s)
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
	implicit real*8 (a-h,o-z)
	dimension T(3,3),Xe(3,1),Xg(3,1),T1(3,3),T2(3,3),
     $	T3(3,3),A(3,3),B(3,3),Ve(3,1),Vg(3,1),P(3,3)
	character*20 inpfile
	data T1,T2,T3 / 9*0d0,9*0d0,9*0d0 /
c
	pi=4d0*datan(1d0)
c
c------Transformation matrix from equatorial to galactic system 
c       (for equatorial coordinates in J2000 system !)
	agp=12.857306*pi/12d0
	dgp=27.128333*pi/180d0
	the=122.93194*pi/180d0
	T1(1,1)=dcos(agp)
	T1(1,2)=dsin(agp)
	T1(2,1)=dsin(agp)
	T1(2,2)=-dcos(agp)
	T1(3,3)=1d0
	T2(1,1)=-dsin(dgp)
	T2(1,3)=dcos(dgp)
	T2(3,1)=dcos(dgp)
	T2(3,3)=dsin(dgp)
	T2(2,2)=-1d0
	T3(1,1)=dcos(the)
	T3(1,2)=dsin(the)
	T3(2,1)=dsin(the)
	T3(2,2)=-dcos(the)
	T3(3,3)=1d0
c
	call Mtxmlt(3,3,3,T3,T2,P)
	call Mtxmlt(3,3,3,P,T1,T)
c
c------Distance of Sun from galactic center in kpc, der IAU Standard:
	Xs=-8.5
c	Xs=-8.0
c
c------Solar motion and motion of LSR in km/s:
c	Vsx=11.0
c	Vsy=5.3
c	Vsz=7.0
	Vlsr=220.
c
c       Altmann Werte
!        Vsx=10.0
!        Vsy=15.0
!        Vsz=7.5
!        Vlsr=225.0

c       Dehnen +Binney 
        Vsx=10.0
        Vsy=5.25
        Vsz=7.17
        Vlsr=220.0
c------Einlesen von Kommandozeile	
c       
c------Reading position, distance, radial velocity and proper motions
c      from the input file (mind the sign of the declination !!)
c
 100	read(*,*) iah,iam,as,idg,idm,ds,dist,rv,pma,pmd
c
	alp=(as/60d0+iam)/60d0 + iah
	alp=alp*pi/12d0
	isig=1
	if(idg.lt.0) then
	 isig=-1
	 idg=-idg
	end if
	del=(ds/60d0+idm)/60d0 + idg
	del=isig*del*pi/180d0
c
	ca=dcos(alp)
	cd=dcos(del)
	sa=dsin(alp)
	sd=dsin(del)
c
	A(1,1)=ca*cd
	A(1,2)=-sa
	A(1,3)=-ca*sd
	A(2,1)=sa*cd
	A(2,2)=ca
	A(2,3)=-sa*sd
	A(3,1)=sd
	A(3,2)=0d0
	A(3,3)=cd
c
	call Mtxmlt(3,3,3,T,A,B)
c
c------Transformation of position vector
c
	do i=1,3
	Xe(i,1)=A(i,1)
	end do
	call Mtxmlt(3,3,1,T,Xe,Xg)
	X=dist*Xg(1,1)+Xs
	Y=dist*Xg(2,1)
	Z=dist*Xg(3,1)
c
c------Transformation of velocity vector
c
	Ve(1,1)=rv
	Ve(2,1)=dist*pma*4.7404
	Ve(3,1)=dist*pmd*4.7404
	call Mtxmlt(3,3,1,B,Ve,Vg)
	U=Vg(1,1)+Vsx
	V=Vg(2,1)+Vsy+Vlsr
	W=Vg(3,1)+Vsz
c
c------Output of galactocentric components in Standardout
c
c	write(2,801) x,y,z
c        write(2,802) U,V,W
        write (*,801) x,y,z,U,V,W
c
 777	format(i2,i6,7x,i2,i3,f5.1,1X,a1,i2,i3,i3,1x,
     $   f4.1,1x,f6.1,13x,f7.2,f7.2)
 801	format(3(f7.3,1X),3(f8.2,1X))
c
	stop 
	end
c
c
c
	SUBROUTINE Mtxmlt(m,l,n,A,B,C)
c	
	implicit real*8 (a-h,o-z)
	dimension A(m,l),B(l,m),C(m,n)
c	
	do i=1,m
	do j=1,n
	C(i,j)=0d0
	do k=1,l
	C(i,j)=C(i,j)+A(i,k)*B(k,j)
	end do
	end do
	end do
	return
	end


