; read pulsation link file and compute and plot nabla and nabla_add
; take care about the centering of these quantities
; note that the lowermost grid point refers to the outer boundary of
; the star.

PRO abla

@physconst

openr,unit,'mod3001.pul',/GET_LUN

i=1l 
readf,unit,i 
a=dblarr(i,32) 
readf,unit,a 
close,unit
free_lun,unit

m=a(1:i-1,1)/XMSUN
na=sqrt(a(0:i-2,21)*a(1:i-1,21))   
n=a(0:i-2,20)

IF 0 THEN BEGIN
    x=n-na
    plot,m,x,/yl,TITLE='20 M_sun star E20@3000',XTITLE='m/M_SUN',YTITLE='abs(N - N_ad)',yr=[1.D-8,1]
    x=na-n
    oplot,m,x,line=1
    
    LEGENDE,/LINKS,/OBEN
    LEGENDE_LINE,0,'N > N_ad'
    LEGENDE_LINE,1,'N < N_ad'
ENDIF ELSE BEGIN
    plot,m,na,TITLE='20 M_sun star E20@3000',XTITLE='m/M_SUN',YTITLE='Nabla'
    oplot,m,n,line=1
    LEGENDE,/LINKS,/UNTEN
    LEGENDE_LINE,0,'N_ad'
    LEGENDE_LINE,1,'N'
ENDELSE

END