Code construction – 1: coordinate grid points and relating parts.

1. Code construction – 1: coordinate grid points and relating parts. Generate coordinate grid points for the spherical coordinate (r, q, f), and tables for functions that is used in the integration of Green’s formula.

2. subroutine coordinate_patch_kit Radial coordinate grid is non-equidistant, angular grids are equidistant. We use the mid-point rule for the quadrature formula. Radial coordinate input parameters: nrg, nrgin, Rin, Rmid, Rout. Rin < Rmid < Rout : r 2 [Rin, Rout], the grid spacing rule is changed at Rmid. nrg : Total number of radial grid points (-1). nrgmid : Number of radial grid points from Rin to Rmid (-1). rg(ir), ir = 0, nrg : Radial coordinate grid points ri . hrg(ir) : Mid-point of radial grids. drg(ir) : Radial grid spacing Dri := ri – ri-1. Dri hri-1hri r0 = Rin = 0 rnrgmid = Rmid » 2R0 rnrg = Rout »104R0 ri-1ri Grid spacing : r 2 [Rin, Rmid], equidistant or let Dr1» 1/4Drnrgmid etc. r 2 [Rmid, Rout], Dri+1 = kDri .

3. Radial Green’s function nrg, rg(ir), hrg(ir), lgmax are used. lgmax: a number of the largest multipole Since mid-point rule is used for the numerical integration, r’ is hrg(irr). hgfn_nb(irr,il,ir) = Dri hri-1hri r0 = Rin = 0 rnrgmid = Rmid » 2R0 rnrg = Rout »104R0 ri-1ri

4. q coordinate input parameter: ntg. q2 [0, p], the grid spacing is equidistant. ntg : Total number of q grid points (-1). thg(it), it = 0, ntg : q grid points qi . hthg(it) : Mid-point of q gridsqi+1/2. dthg : q grid spacing Dq. Functions associated with q coordinate input parameters: lgmax. Trigonometric funcitons. sinthg(it), it = 0, ntg : sinqi. hsinthg(it) : sinqi+1/2 sin at mid-point of q grids. costhg(it), it = 0, ntg : cosqi. hcosthg(it) : cosqi+1/2 cosin at mid-point of q grids. Associated Legendre funcitons plmg(it,il,im), it = 0, ntg , il = 0, lgmax , im = 0, lgmax : at each grid point qi hplmg(it,il,im), it = 1, ntg , il = 0, lgmax , im = 0, lgmax : at each half grid point qi+1/2 Since mid-point rule is used for the numerical integration, we need hplmg(it,il,im) at half grid points qi+1/2.

5. f coordinate input parameter: npg. f2 [0, 2p], the grid spacing is equidistant. npg : Total number of f grid points (-1). phig(ip), ip = 0, npg : f grid points fi . hphig(ip) : Mid-point of f gridsfi+1/2. dphig : f grid spacing Df. Functions associated with f coordinate input parameters: mgmax = lgmax. Trigonometric funcitons. sinphig(ip), ip = 0, npg : sinfi. hsinphig(ip) : sinfi+1/2 sin at mid-point of f grids. cosphig(ip), ip = 0, npg : cosfi. hcosphig(ip): cosfi+1/2 cos at mid-point of f grids. Trigonometric funcitons with the angle mf. sinmpg(im,ip), im = 0, mgmax , ip = 0, npg : at each grid point fi cosmpg(im,ip), im = 0, mgmax , ip = 0, npg : at each grid point fi hsinmpg(im,ip), im = 0, mgmax , ip = 1, npg : at each half grid point fi+1/2 hcosmpg(im,ip), im = 0, mgmax , ip = 1, npg : at each half grid point fi+1/2

6. Weight for the integration: assigned at the mid-points. wrg(irr) = hrg(irr)2 drg(irr) : weight for the radial integratoin ri2Dri wtg(itt) = hsintheg(itt) dtg : weight for the q integratoin sinqi Dqi wpg(ipp) = dpg : weight for the f integratoin Dfi Weights for the numerical integration formula of mid-point rule is included in the above (which is 1 at all points). Other parts in the Green’s formulainput parameters: lgmax = mgmax. faclmg(il,im) = il = 0, lgmax, im = 0, mgmax. epsig(im) = end subroutine coordinate_patch_kit