;
; ghv_csss:  Calculate harmonic coefficients gcp & hcp using the calculated
;            vector magnetic field at CP surface and the CSSS model
; input:     apar,Nmax,Rcp,Rss,aith,ajph,sth,cm,sm,P,dP,bbr,bbt,bbp,ghcpf
;
; output:    ghcpf,gcp,hcp
; 
; reproduced: Mar. 5, 96
; MODIFIED:   Mar. 13, 2006	Xuepu Zhao for Nagoya group
;
pro ghv_csss,bbr,bbt,bbp,sth,cm,sm,P,dP,apar,Rcp,Rss,gcp,hcp,ghcpf

    if N_params() LT 1 then begin
       print,'ghv_csss,bbr,bbt,bbp,sth,cm,sm,P,dP,apar,Rcp,Rss,gcp,hcp,ghcpf'
       print,'INPUT: bbr,bbt,bbp,sth,cm,sm,P,dP,apar,Rcp,Rss'
       print,'OUTPUT: gcp,hcp,ghcpf
       return
    endif
print,'This is /home/xuepu/ZIDLLIB/CSSS/ghv_csss.pro'    
    sz=SIZE(bbr) & iph=sz(1) & jth=sz(2)
    sz=SIZE(cm) & nmax=sz(1)-1
; Determine numbers of elements for various demensions of arrays
    Iph = fix(iph + 0.5)
    Jth = fix(jth + 0.5)
    Ixx = (nmax+1)^2
    Iyy = fix(3. * Iph * Jth + 0.5)
; Define arrays
    Bout = fltarr(Iyy)
    ab = fltarr(Iyy,Ixx)
    GHout = fltarr(Ixx)
    gcp = fltarr(nmax+1,nmax+1)
    hcp = gcp
; calculate ab_ixiy (atrray ab(Iyy,Ixx)). iy--i,j,k; ix--n,m
    ar1=fltarr(nmax+1)
    ar2=ar1
    ar3=fltarr(nmax+1,nmax+1)
    rssa = Rss+apar
    rca = Rcp+apar
    FOR n=0, nmax do begin
        l = 2*n+1
        rssal = rssa^l
        rcal = rca^l
        Kn = (rssal-rcal)*Rcp/rca
        Kn = Kn/((n+1)*rssal + n*rcal)
        ar1(n) = 1.0
        ar2(n) = -1.*Kn
        FOR m=0, n do ar3(n,m) = m*Kn
    ENDFOR
    ix = -1
    FOR n=0, nmax do begin
        FOR m=0, n do begin
            ix = ix + 1
            iy = 0
            FOR j=0, Jth-1 do begin
                FOR i=0, Iph-1 do begin
                    ab(iy,ix) = ar1(n) * cm(m,i) * P(n,m,j)
                    iy = iy + 1
                    ab(iy,ix) = ar2(n) * cm(m,i) * dP(n,m,j)
                    iy = iy + 1
                    ab(iy,ix) = ar3(n,m) * sm(m,i) * P(n,m,j) / sth(j)
                    iy = iy + 1
                ENDFOR
            ENDFOR
        ENDFOR
    ENDFOR
; print,'ix=',ix
    FOR n=1, nmax do begin
        FOR m=1, n do begin
            ix = ix + 1
            iy = 0
            FOR j=0, Jth-1 do begin
                FOR i=0, Iph-1 do begin
                    ab(iy,ix) = ar1(n) * sm(m,i) * P(n,m,j)
                    iy = iy + 1
                    ab(iy,ix) = ar2(n) * sm(m,i) * dP(n,m,j)
                    iy = iy + 1
                    ab(iy,ix) = -ar3(n,m) * cm(m,i) * P(n,m,j) / sth(j)
                    iy = iy + 1
                ENDFOR
            ENDFOR
        ENDFOR
    ENDFOR
; get abt=ab^{-1}=abt(Ixx,Iyy), abba=array(Ixx,Ixx)
    abt = TRANSPOSE(ab)
    abba = abt # ab
    abbai = INVERT(abba)
   ;help,ab,abt,abba,abbai
; calculate Bout(Iyy). iy -- i,j,k (theta, phi, rtp)
    Bonly = 0
    ss = 0
    iy = 0
    FOR j=0, Jth-1 DO BEGIN
        FOR i=0, Iph-1 do begin
            if (bbr(i,j) lt 0.) then begin
               Bout(iy) = -1 * bbr(i,j)
               iy = iy + 1
               Bout(iy) = -1 * bbt(i,j)
               iy = iy + 1
               Bout(iy) = -1 * bbp(i,j)
; For checking the program with pure potential-like field
;           if (bbr(i,j) lt 0.) then begin
;              Bout(iy) = bbr(i,j)
;              iy = iy + 1
;              Bout(iy) = bbt(i,j)
;              iy = iy + 1
;              Bout(iy) = bbp(i,j)
            endif else begin
               Bout(iy) = bbr(i,j)
               iy = iy + 1
               Bout(iy) = bbt(i,j)
               iy = iy + 1
               Bout(iy) = bbp(i,j)
            endelse
            iy = iy + 1
        ENDFOR
    ENDFOR
    abB = Bout # ab
    GHout = abB # abbai
; get gcp and hcp and create a data file for g & h at Rcp
    ix = -1
    FOR n=0, nmax do FOR m=0, n do begin
                         ix = ix + 1
                         gcp(n,m) = GHout(ix)
                     ENDFOR
    FOR n=1, nmax do FOR m=1, n do begin
                         ix = ix + 1
                         hcp(n,m) = GHout(ix)
                     ENDFOR
    openw,10,ghcpf
    printf,10
    printf,10,nmax
    FOR n=0, nmax do FOR m=0, n do printf,10,n,m,gcp(n,m),hcp(n,m)
    close,10

end