ngc	# Rebecca Stanek
	# ASTR 327
	expand 1.24
	ctype default

	WINDOW 2 2 1 2 
	load math
	data ngc2403.dat
	read { r 2 vc 3 err 4 }
	set r = 1000*r
	set vcsq = vc*vc
	define h 6800
	set y = (r/(2*$h))
	define G 4.301E-3
	define sig 18.4
	lim r vcsq
	box
	connect r vcsq
	#set mess = (-.38407-ln(y))*y**2 + (.399449-.75*ln(y))*y**4 + (.156135-.15625*ln(y))*y**6
	set mess = y*y*(I0(y)*K0(y)-I1(y)*K1(y))
	set model = (4*PI*$G*$sig*$h)*(mess)
	ctype yellow
	connect r model
	ctype default
	xlabel R
	ylabel V_c^2

	WINDOW 2 2 2 2 
	set DIMEN(diff) = 33
	do i = 0, 32 {set diff[$i] = vcsq[$i] - model[$i]}
	#lim r 0 10E7 
	#box
	#connect r diff
	
	set mass = (diff*r)/$G
	lim r mass
	connect r mass
	ctype default
	box
	ctype default
	xlabel R (pc)
	ylabel M(R) (Solar Masses)

	WINDOW 2 2 1 1 
	lsq r mass 
	set rho = (1/(4*PI*r*r))*$a
	lim r 0 .75 
	box
	connect r rho
	xlabel R
	ylabel \rho(R)

	WINDOW 2 2 2 1
	set logrho = ln(rho)
	set logr = ln(r)
	lim logr logrho 
	box
	connect logr logrho
	set DIMEN(x2) = 33
	set DIMEN(y2) = 33
	lsq logr logrho x2 y2 
	#lim x2 y2
	echo $a 
	set xfit = 2, 12
	set yfit = $a*xfit + $b
	ctype cyan
	connect xfit yfit
	ctype default	
	xlabel ln(R)
	ylabel ln(\rho)
	relocate 7 3
	label \alpha = 2