if !d.name eq 'PS' then begin
  device,xsize=18,ysize=9,yoffset=3
  !p.charthick=2.4 & !p.thick=2.4 & !x.thick=2.4 & !y.thick=2.4
end
;
;  mv idl.ps ../fig/pdet.ps
;  mv idl.ps ../fig/pdet_large.ps
;
!p.charsize=1.4
!x.margin=[7.5,0.6]
!y.margin=[3.2,0.4]
!p.multi=0
!x.title='!8k!d!9x!n!6'
!y.title='!6Ra'
circ_sym,1.3,1
;
ii=complex(0.,1.)
onethird=1./3.
sq3=sqrt(3.)
Rac=1707.761d0
kpc=3.117d0
;
a1=-1.
a2=+.5+.5*ii*sqrt(3.)
a3=+.5-.5*ii*sqrt(3.)
;
;  compute the 3 values of q
;
nRa=201 & nkp=301
;
RRa=grange(100.,3000.,nRa) & kkp=grange(.5,5.0,nkp) & levx=50
RRa=grange(1680.,1800.,nRa) & kkp=grange(2.5,3.8,nkp) & levx=1.
;
det=fltarr(nkp,nRa)
;
for iRa=0,nRa-1 do begin
for ikp=0,nkp-1 do begin
  kp=kkp[ikp]
  Ra=RRa[iRa]
  ;
  kp4=kp^4
  ;
  q0=kp*sqrt((Ra/kp4)^onethird-1.)
  q1=kp*sqrt(1.+a1*complex(Ra/kp4)^onethird)
  q2=kp*sqrt(1.+a2*complex(Ra/kp4)^onethird)
  q3=kp*sqrt(1.+a3*complex(Ra/kp4)^onethird)
  ;
  c1=(q1^2-kp^2)^2
  c2=(q2^2-kp^2)^2
  c3=(q3^2-kp^2)^2
  ;
  ;m=[[   cosh(.5*q1),   cosh(.5*q2),   cosh(.5*q3)], $
  ;   [q1*sinh(.5*q1),q2*sinh(.5*q2),q3*sinh(.5*q3)], $
  ;   [c1*cosh(.5*q1),c2*cosh(.5*q2),c3*cosh(.5*q3)] ]
  ;
  ;det(ikp,iRa)=LA_determ(m)
  ;det=imaginary((sq3+ii)*q2*tanh(.5*q2))+q0*tan(.5*q0)
  det(ikp,iRa)=imaginary((sq3+ii)*q2*tanh(.5*q2))+q0*tan(.5*q0)
endfor
endfor
;
lev0=[-1,0,1]*1e-6
lev=grange(-1,1,40)*levx
print,minmax(det)
;lev=grange(min(abs((det))),max(abs((det))),20)
contour,clip(det,minmax(lev)),kkp,RRa,lev=lev,/fil
contour,det,kkp,RRa,lev=lev0,/over,col=255
oplot,[1,1]*kpc,[1,1]*Rac,col=55,ps=8
;print,kp,Ra,abs(det)
END