      parameter(idim=6500000)
      character*400 filename
      dimension rm(idim),x(idim),y(idim),z(idim),vx(idim),vy(idim),vz(idim)
      dimension pot(idim)
       dimension rx(idim),ry(idim),rz(idim)
       dimension xn(idim),yn(idim),zn(idim)
       dimension vxn(idim),vyn(idim),vzn(idim)


      nstr=1000000

      write(filename,'(''../data/satellite.dat'')')
      write(filename,'(''../data/isolation.dat'')')
      open(1,file=filename)
      do i=1,nstr
      read(1,*)rm(i),rx(i),ry(i),rz(i),vx(i),vy(i),vz(i),pot(i)
      enddo
	ngas=0
	ndrk=0

	munit=2.3262e5
	rcut=1e10
	npar=ngas+ndrk+nstr
      write(filename,'(''../data/cma.dat'')')
	open(1,file=filename)
        read(1,1)xcm,ycm,zcm

      write(filename,'(''../data/vcma.dat'')')
	open(1,file=filename)
        read(1,*)vxcm,vycm,vzcm

	rjx=0.
	rjy=0.
	rjz=0.
	do j=1,nstr
	i=j+ngas+ndrk
      rxi=rx(i)-xcm
      ryi=ry(i)-ycm
      rzi=rz(i)-zcm
      vxi=vx(i)-vxcm
      vyi=vy(i)-vycm
      vzi=vz(i)-vzcm
	r=sqrt(rx(i)**2+ry(i)**2+rz(i)**2)
	if(r.le.rcut)then
	rjx=rjx+rm(i)*(ryi*vzi-rzi*vyi)
	rjy=rjy+rm(i)*(rzi*vxi-rxi*vzi)
	rjz=rjz+rm(i)*(rxi*vyi-ryi*vxi)
	endif
	enddo
	close(1)
	write(*,*)rjx/1e10,rjy/1e10,rjz/1e10

	rj=sqrt(rjx**2+rjy**2+rjz**2)
	rjx=rjx/rj
	rjy=rjy/rj
	rjz=rjz/rj
	e1x= rjy*rjz
	e1y=-rjx*rjz
	e1z=0.
	e1=sqrt(e1x**2+e1y**2+e1z**2)
	e1x=e1x/e1
	e1y=e1y/e1
	e1z=e1z/e1

	e2x=  rjx*rjz**2
	e2y=  rjy*rjz**2
	e2z=-(rjx**2+rjy**2)*rjz
	e2=sqrt(e2x**2+e2y**2+e2z**2)
	e2x=e2x/e2
	e2y=e2y/e2
	e2z=e2z/e2

	e3x=rjx
	e3y=rjy
	e3z=rjz
	e3=sqrt(e3x**2+e3y**2+e3z**2)
	e3x=e3x/e3
	e3y=e3y/e3
	e3z=e3z/e3

      write(filename,'(''../data/exyza.dat'')')
        open(1,file=filename)
	write(1,*)e1x,e1y,e1z
	write(1,*)e2x,e2y,e2z
	write(1,*)e3x,e3y,e3z
	close(1)


 	do i=1,npar
       rx(i)=rx(i)-xcm
       ry(i)=ry(i)-ycm
       rz(i)=rz(i)-zcm
       vx(i)=vx(i)-vxcm
       vy(i)=vy(i)-vycm
       vz(i)=vz(i)-vzcm
	rm(i)=rm(i)*munit
	enddo

c       realiza la rotacion de las posiciones  velocidades xnew=a_transpuesta*xold
	do i=1,npar
	xn(i)=e1x*rx(i)+e1y*ry(i)+e1z*rz(i)
	yn(i)=e2x*rx(i)+e2y*ry(i)+e2z*rz(i)
	zn(i)=e3x*rx(i)+e3y*ry(i)+e3z*rz(i)
	vxn(i)=e1x*vx(i)+e1y*vy(i)+e1z*vz(i)
	vyn(i)=e2x*vx(i)+e2y*vy(i)+e2z*vz(i)
	vzn(i)=e3x*vx(i)+e3y*vy(i)+e3z*vz(i)
	enddo

c       graba masas y posiciones gas
      write(filename,'(''../data/rota200_gas.dat'')')
	open(1,file=filename)
	do j=1,ngas
	i=j
      write(1,1)rm(i),xn(i),yn(i),zn(i),vxn(i),vyn(i),vzn(i),-pot(i)
	enddo
	close(1)

	write(*,*)ndrk
c       graba masas y posiciones dark 
      write(filename,'(''../data/rota200_drk.dat'')')
	open(1,file=filename)
	do j=1,ndrk
	i=ngas+j
      write(1,1)rm(i),xn(i),yn(i),zn(i),vxn(i),vyn(i),vzn(i),-pot(i)
	enddo
	close(1)

c       graba masas y posiciones estrellas
      write(filename,'(''../data/rota200_str.dat'')')
	open(1,file=filename)
	do j=1,nstr
	i=ngas+ndrk+j
      write(1,1)rm(i),xn(i),yn(i),zn(i),vxn(i),vyn(i),vzn(i),-pot(i)
	enddo
	close(1)

 1      format(30e18.8)
	end
