c
	program nstars_sc
c
c	NSTARS program with the sun in the center in the output
c
c-------------------------------------------------------------------
c
	implicit none
c
	include 'galactic.fi'
c
c----- input variables
c
	character name1*9
c
	real ras
	integer rah, ram, decd, decm, decs
	character decsign*1
c
	real parallax
c
c----- calculation variables
c
	real ra, dec, e_ra, e_dec
	real epoch
	real ra_std, dec_std
c
	real lii, bii, e_lii, e_bii
c
	real distance_pc, e_distance_pc
c
	real x, y, z, e_x, e_y, e_z
c
	character input_file*30, output_file*30
c
	real xcenter, ycenter, zcenter
	real dx, dy, dz
	real r, cut_radius_pc
c
c-------------------------------------------------------------------
c
c----- Cut radius set to 5 parsec diameter
c
	xcenter = sun_galactic_x_pc
	ycenter = sun_galactic_y_pc
	zcenter = sun_galactic_z_pc
c
	cut_radius_pc = 5.0
c
	input_file = 'nstars.dat'
	output_file = 'nstars_sc.xyz'
c
	open(11,file=input_file,status='old')
	open(12,file=output_file,status='unknown')
	rewind(12)
	open(13,file='nstars_sc_radii.dat',status='unknown')
	rewind(13)
c
	write(12,1100)
 1100	format(/'       Relative Positions of Stars in the Solar Neighborhood'//
     &       ' -- positions scaled from -1 to 1'/
     &       ' -- all stars fit within sphere of unit radius'/
     &       ' -- true sphere scale is 5 parsecs in radius (10 parsec diameter)'///
     &       ' Name               X          Y          Z'/
     &       ' -----------      -----      -----      -----')
c
c
	read(11,1009)
 1009	format(/)
c
  100	read(11,1010,end=900) name1, rah, ram, ras, decsign, decd,
     &			      decm, decs, parallax
 1010	format(a9,43x,i2,1x,i2,1x,f4.1,6x,a1,i2,1x,i2,1x,i2,39x,f7.5)
c
	print*,name1
c
c----- Calculate RA in decimal hours and DEC in decimal degrees
c
	ra = float(rah) + float(ram)/60.0 + ras/3600.0
c
	dec = float(decd) + float(decm)/60.0 + float(decs)/3600.0
	if (decsign .eq. '-') dec = -1.0*dec
c
c----- Calculate distance in parsecs
c
	if (parallax .gt. 0.0) then
	  distance_pc = 1.0/parallax
	else
	  print*,' Non-positive parallax - ',name1,' - ',parallax
	  distance_pc = 1.0/206265.0
	endif
c
c----- No error bars given for this data
c
	e_ra = 0.0
	e_dec = 0.0
	e_distance_pc = 0.0
c
c----- Assume epoch is J2000.0
c
	epoch = 2000.0
c
c----- Convert from celestial coordinates to galactic coordinates
c
	call celestial_to_galactic(ra, dec, e_ra, e_dec, epoch, ra_std,
     &				   dec_std, lii, bii, e_lii, e_bii)
c
c----- Convert from galactic coordinates and distance to 3D coordinates
c
	call galactic_to_xyz(lii, bii, e_lii, e_bii, distance_pc,
     &			     e_distance_pc, x, y, z, e_x, e_y, e_z)
c
c----- Calculate distance from center
c
	dx = x - xcenter
	dy = y - ycenter
	dz = z - zcenter
c
	r = sqrt(dx*dx + dy*dy + dz*dz)
c
c----- Scale and shift coordinates
c
ccc	dx = (dx + cut_radius_pc) / (2.0 * cut_radius_pc)
ccc	dy = (dy + cut_radius_pc) / (2.0 * cut_radius_pc)
ccc	dz = (dz + cut_radius_pc) / (2.0 * cut_radius_pc)
c
	dx = dx/cut_radius_pc
	dy = dy/cut_radius_pc
	dz = dz/cut_radius_pc
c
c----- Write out coordinates
c
	if (r .le. cut_radius_pc) then
	  write(12,1020) name1, dx, dy, dz
 1020	  format(1x,a11,3(1x,f10.3))
c
	else
	  print*,' Star outside cut radius - ',name1
	  print*,'   radius = ',r
	endif
c
	write(13,1030) name1, r
 1030	format(1x,a11,1x,f11.4)
c
	goto 100
c
  900	continue
c
	close(11)
	close(12)
	close(13)
c
	end