function vr, l, palm, r
bigg = 4.43e-6
vr = ((!pi*bigg*palm*l)/(2*r))
return, vr
end

function lum, r
rs = 6.25
sigma0 = 80.7e6
l = sigma0 * 2 * !pi * rs * (rs - exp(-r/rs)*(r+rs))
return, l
end

function massr, mnot, cr, r
mr = 4*!pi*mnot*((-r/(r+cr)) + alog((cr+r)/cr))
return, mr
end

function vhalo, r, mr
bigg = 4.43e-6
vh = bigg * mr / r
return, vh
end

pro gravity, palm, cm, cr

 device, decomposed=0
  red=  [1,0,1,0,0,1,1,0]
  green=[1,0,0,1,0,0,1,1]
  blue= [1,0,0,0,1,1,0,1]
  tvlct,255*red,255*green,255*blue

mnot = cm/(4*!pi*(-0.5+alog(2.)))
r = [0.00001, .325, .974, 1.95, 2.6, 3.57, 7.47, 12.3, 17.2, 22.1, 27.0, 31.8, 36.7, 41.6, 46.5]
lvect = findgen(15)
vrvect = findgen(15)
mrvect = findgen(15)
vnvect = findgen(15)
l = lum(r)
mrvect = massr(mnot, cr,r)
print, mrvect
vhvect = vhalo(r, mrvect)
print, vhvect

vreal = [0.0, 90, 125, 148, 180, 201, 206, 214, 212, 213, 214, 213, 213, 203, 198]

for i = 0, 14 do begin
        l = lum(r[i])
        lvect[i] = l
        vrvect[i] = vr(l, palm, r[i])
endfor
vsvect = sqrt(vrvect)
vnvect = sqrt (vrvect + vhvect)
window, 0, title='palm tree plot thing'
        plot,r,vsvect, ytitle = 'rotational Velocity (km/s)', xtitle = 'Radius (kpc)', background=0, color=1, psym=2
        oplot, r, vreal, psym=5, color=3
        oplot, r, vnvect, psym = 5, color = 5
end