pro lmbdrk_sbsp,scnpos,srad,xctr,yctr,lmdrk

;  routine to compute the theoretical limb darkening function for
;  prescribed slit scan positions, solar radius, and center of the
;  sp scan map

;  INPUTS:	scnpos	= linear array of the scan center positions
;		srad	= solar radius at time of observation, arcsec
;		xctr	= x coordinate of map center, arcsec
;		yctr	= y coordinate of map center, arcsec
;  OUTPUTS:	lmdrk	= 2-D array of intensity of limb darkening 
;			  relative to disk center for prescribed slit
;			  scan positions

;  This assumes these are normal maps
;  Limb darkening function from Neckel and Labs, SP, 153, p.97, 1994

;  limb darkening polynomial coefficients vs mu
;  values of aa are for wavelength 640.97 nm
aa = [0.33644, 1.30590, -1.79238, 2.45040, -1.89979, 0.59943]

nnx = n_elements(scnpos)

;  generate arrays of map coordinates in arcseconds
;  use the nominal values for step size
;xxpos = -1023. + findgen(2047)
xxpos = scnpos*0.1476
yypos = -511. + findgen(1024)
yypos = yypos*.1585
ypos = replicate(1.,nnx)#yypos + yctr
xpos = xxpos#replicate(1.,1024) + xctr
xpos = double(xpos)
ypos = double(ypos)


npoly = n_elements(aa)
;srad = 970.	; try sample solar radius in arcsec
;  find arcsec radius to each point in map
radd = sqrt(xpos*xpos + ypos*ypos)
radf = radd/srad
mu = sqrt(1.-radf*radf)
lmdrk = mu
lmdrk(*,*) = 0.D
for kk=0,npoly-1 do lmdrk = lmdrk+aa(kk)*mu^kk
;tvwinp,lmdrk

;stop
end