;
; For HCCS and HCCSSS calculation below r=Rcp
; H(n,..)={(1+a)^n/(n+1)}/(a+r)^{(n+1)}  n=0,1,...Nmax
; dHr(n,,..)=dH(n,..)/dr=-{(1+a)^n/(n+1)}{(n+1)/(r^2*(r+a)^n)}
;
pro hchf,Nmax,apar,r,Hr,dHr,d2Hr

    Hr = fltarr(Nmax+1)
    dHr = fltarr(Nmax+1)
    a1 = 1.0 + apar
    ar = r + apar
    for n = 0,Nmax do begin
        n1 = n + 1
        un = a1^n/n1
        Hr(n) = un / ar^n1
        dHr(n) = -un * n1 /(r^2.*ar^n)   ; Modified 10 March 1998
    endfor

end