;
; ghv_csss:  Calculate g & h using a regular synoptic map of Vector
;             magnetic field and the CSSS model
; input:     apar,Nmax,Rcp,Rss,aith,ajph,sth,cm,sm,P,dP,bbr,bbt,bbp,ghcsssf
;
; output:    ghcsssf,gcp,hcp
;
; reproduced: Mar. 5, 96
;
pro ghv_csss,apar,Nmax,Rcp,Rss,aith,ajph,sth,cm,sm,P,dP,bbr,bbt,bbp,$
             ghcsssf,gcp,hcp

;tjd+
	print, 'computing...'
	print, 'computing...'
	print, 'go get a cup of coffee...'
	print, 'computing...'
;tjd-

    if N_params() LT 1 then begin
       print,'SYNTAX - ghv_csss,apar,Nmax,Rcp,Rss,aith,ajph,sth,cm,sm,'
       print,'                       P,dP,bbr,bbt,bbp,ghcsssf,gcp,hcp'
       return
    endif
; Determine numbers of elements for various demensions of arrays
    Ith = fix(aith + 1 + 0.5)
    Jph = fix(ajph + 1 + 0.5)
    Iyy = fix(3. * Jph * Ith + 0.5)
  ; Iyy = 3. * Jph * Ith
    Ixx = (Nmax+1)^2
; 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 i=0, Ith-1 do begin
                FOR j=0, Jph-1 do begin
                    ab(iy,ix) = ar1(n) * cm(m,j) * P(n,m,i)
                    iy = iy + 1
                    ab(iy,ix) = ar2(n) * cm(m,j) * dP(n,m,i)
                    iy = iy + 1
                    ab(iy,ix) = ar3(n,m) * sm(m,j) * P(n,m,i) / sth(i)
                    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 i=0, Ith-1 do begin
                FOR j=0, Jph-1 do begin
                    ab(iy,ix) = ar1(n) * sm(m,j) * P(n,m,i)
                    iy = iy + 1
                    ab(iy,ix) = ar2(n) * sm(m,j) * dP(n,m,i)
                    iy = iy + 1
                    ab(iy,ix) = -ar3(n,m) * cm(m,j) * P(n,m,i) / sth(i)
                    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 i=0, Ith-1 DO BEGIN
        FOR j=0, Jph-1 do begin
            if (bbr(j,i) lt 0.) then begin
               Bout(iy) = -1 * bbr(j,i)
               iy = iy + 1
               Bout(iy) = -1 * bbt(j,i)
               iy = iy + 1
               Bout(iy) = -1 * bbp(j,i)
; For checking the program with pure potential-like field
;           if (bbr(j,i) lt 0.) then begin
;              Bout(iy) = bbr(j,i)
;              iy = iy + 1
;              Bout(iy) = bbt(j,i)
;              iy = iy + 1
;              Bout(iy) = bbp(j,i)
            endif else begin
               Bout(iy) = bbr(j,i)
               iy = iy + 1
               Bout(iy) = bbt(j,i)
               iy = iy + 1
               Bout(iy) = bbp(j,i)
            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
  ;  print,'the mystery ',ghcsssf ;tjd added print
    openw,10,ghcsssf
    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