;
; For HCSS calculation below r=Rss
; U=[1/(a+1)]^n*{n+1+n[(a+1)/(a+Rss)]^(2n+1)}
; H(n,..)=[1/(a+r)]^(n+1)*{1-[(a+r)/(a+Rss0]^(2n+1)}/U
;        =[1/(a+r)]*[(a+1)/(a+r)]^n*{1-[(a+r)/(a+Rss)]^(2n+1)}/{n+1+n[(a+1)/(a+Rss)]^(2n+1)}
; dH(n,..)/dr=-[1/r]^2*[1/(r+a)]^n*{n+1+n*[(a+r)/(a+Rss)]^(2n+1)}/U
;        =-[1/r]^2*[(a+1)/(r+a)]^n*{n+1+n*[(a+r)/(a+Rss)]^(2n+1)}/{n+1+n[(a+1)/(a+Rss)]^(2n+1)}
; ar=a+r; a1arn={(a+1)/(a+r)}^n; ararsn={(a+r)/(a+Rss)}^{2n+1}
; a1arsn=n+1+n*{(a+1)/(a+Rss)}^{2n+1}
; H(n,..)={a1arn/ar}{1-ararsn}/a1arsn
; dHr(n,,..)=dH(n,..)/dr=-{a1arn/r^2}{n+1+n*ararsn}/a1arsn
; d2Hr(n,..)={a1arn/ar^3}{(n+1)(n+2)-n(n-1)ararsn}/a1arsn
; MODIFIED: 	28May2003 	Xuepu Zhao for high Nmax
; MODIFIED:     01APR2004	Xuepu Zhao for adding /dbl
;
pro zsshf,nmax,rss,r,Hr,dHr,d2Hr,DBL=dbl

    if keyword_set(dbl) then Hr = dblarr(nmax+1) $
                        else Hr = fltarr(nmax+1)
    dHr = Hr
    d2Hr = Hr
    apar=0.0
    a1 = 1 + apar
    ar = r + apar
    ars = rss + apar
    a1ar = double(a1/ar)
    a1ars = double(a1/ars)
    arars = double(ar/ars)
    for n = 0,nmax do begin
        n1 = n + 1.0
        n2 = n + 2.0
        n3 = n + 3.0
        n21 = 2*n + 1.0
        a1arn = a1ar^n
        a1arsn = n1 + n * a1ars^n21
        ararsn = arars^n21
        Hr(n) = (a1arn/ar) * (1 - ararsn) / a1arsn
        dHr(n) = -(a1arn/r^2) * (n1+n*ararsn) / a1arsn
    ;   d2Hr(n) = (a1arn/ar^3)*(n1*n2-n*(n-1)*ararsn)/a1arsn
    endfor

end