; this is a code to compute the stability function

pro psimphih,h,mo,psim,phih

beta=5.
gamma=16.0

IF mo GE 0 THEN BEGIN

phih=(1+beta*h/mo)
psim=-beta*h/mo

ENDIF ELSE BEGIN

phih=1.0/sqrt(1-gamma*h/mo)
y=(1-gamma*h/mo)^(0.25)

psim=2*alog((1+y)/2.0)+alog((1+y*y)/2)+2*atan((1-y)/(1+y))

ENDELSE

return
END