pro jam_rms, mbh, ml1, log_dm,rhalo, beta2, inc, rmsmod, XS=xs, YS=ys

if keyword_set(xs) then if xs mod 2 ne 1 then message,$
 'Set Odd Xsize' else xsize=xs else xsize=101
if keyword_set(ys) then if ys mod 2 ne 1 then message,$
 'Set Odd Ysize' else ysize=ys else ysize=11

sol=mrdfits('boschmge_new.fits',0)
dist=104.2
pix=0.247
peak=10.^(log_dm)/(2.*!pi*(rhalo)^2)
scl=fltarr(n_elements(sol[*,0])+1,n_elements(sol[0,*]))
for qq=0,n_elements(sol[*,0])-1 do scl[qq,*]=sol[qq,*]
scl[n_elements(sol[*,0]),*]=[peak,rhalo/(dist/.206265),0.99,0.]
surf_lum=scl[*,0]
sigma_lum=scl[*,1]
qobs_lum=scl[*,2]
surf_pot=surf_lum
sigma_pot=sigma_lum
qobs_pot=qobs_lum
print,peak

x1=pix*(indgen(xsize)-(xsize-1)/2)
y1=pix*(indgen(ysize)-(ysize-1)/2)


x2=fltarr(xsize,ysize)
y2=x2

m_bh=mbh
s1=where(qobs_lum eq min(qobs_lum),complement=s2) 

;mask_pa=168d
;if keyword_set(pa) then alph=(mask_pa-pa)/180*!PI else alph=!PI
alph=!PI

for i=0,xsize-1 do begin
for j=0,ysize-1 do begin
x2[i,j]=x1[i]*cos(alph)-y1[j]*sin(alph)
y2[i,j]=x1[i]*sin(alph)+y1[j]*cos(alph)
endfor
endfor


sigmapsf=1./2.35
;normpsf=[0.5777,0.4223]

mbh=m_bh

surf_pot[s1]=surf_lum[s1]*ml1
surf_pot[s2]=surf_lum[s2]*ml1

beta_=beta2

print,surf_pot/surf_lum

jam_axisymmetric_rms,surf_lum,sigma_lum,qobs_lum,surf_pot,sigma_pot,qobs_pot,inc,mbh,dist,x2,y2,rmsmodel,PIXSIZE=pix,SIGMAPSF=sigmapsf,STEP=pix/8,BETA=beta_;, GOODBINS=good, RMS=rms,ERMS=e_rms


rmsmod=rmsmodel




print,'==========================DONE============================'

end