function dark_factor, sep, param, qstop=qstop
sep = [sep]
if n_elements(param) eq 0 then param = 0.6
if param ge 1 then begin
  print,'Limb-darkening parameter cannot exceed 1, sorry'
  return,1.
endif
n_sep = n_elements(sep)
result = fltarr(n_sep)
img = fltarr(200,200)+1.
x = indgen(200)-100
y = x
xx = fltarr(200,200)
yy = xx
rr = xx
for i = 0,199 do xx(*,i) = x/100.
for i = 0,199 do yy(i,*) = y/100.
disk = where(xx^2+yy^2 lt 1)
rr(disk) = sqrt(xx(disk)^2+yy(disk)^2)
for i = 0, n_sep-1 do begin
  if sep(i) ge 1.99 then result(i) = 1. else begin
  disk2 = where((xx-sep(i))^2+yy^2 lt 1.)
  rr2 = fltarr(200,200)
  rr2(disk2) = sqrt((xx(disk2)-sep(i))^2+yy(disk2)^2)
  mu = fltarr(200,200)
  mu(disk) = sqrt(1-rr(disk)^2)
  dark = fltarr(200,200)
  dark(disk) = param*mu(disk) + (1-param)
  ylim = sep(i)/2
  occ = where(rr2 ne 0 and rr ne 0)
  dark_sum = total(dark(occ))
  dark_norm= total(dark)
  result(i) = 1.-dark_sum/dark_norm
endelse
endfor
if n_elements(qstop) ne 0 then stop
return, result
end