;+ ; NAME : CFF1N (procedure) ; PURPOSE : ; calculate current-free field based on the longitudinal ; field (BL) in rdata (planar geometry) ; CATEGORY : ; solar magnetic field computation package (MAGPACK2) ; CALLING SEQUENCE : ; cff1n, Bl, bxp, byp ; INPUTS : ; Bl -- longitudinal magnetic field ; OPTIONAL INPUT PARAMETERS : ; none ; KEYWORD PARAMETERS : ; position : position on the solar disk, [ ew, ns ] ; disk center => [0,0], limb => +-1. ; dx : scale of 1 pixel in unit of 10000km ; OUTPUTS : ; bxp, byp : current-free transverse fields ; COMMON BLOCKS : none ; SIDE EFFECTS : mparam will be modified ; RESTRICTIONS : none ; PROCEDURE : ; use FFT to calculate current-free fields ; MODIFICATION HISTORY : ; T.Sakurai, 1992 ; K.Ichimoto, 92/10/17 modify cff1.pro ; ;---------------------------------------------------------------------- ; pro setbtp, kwx,kwy, alx,aly, cx,cy, amat ; calculate complex coefficients (cx,cy) to be multiplied ; on the Fourier transform of BL ; kwx,kwy = integer indices ; alx = x-width /2*pi ; amat = matrix of coordinate transform between data and tangential ; coordinats ; ci = complex (0.0, 1.0) akwx=kwx/alx akwy=kwy/aly akx=akwx*amat(0,0) aky=akwx*amat(0,1)+akwy*amat(1,1) akz=sqrt(akx*akx+aky*aky) cc=ci*akz-akwx*amat(0,2)-akwy*amat(1,2) cx=(akwx+cc*amat(0,2))/(cc*amat(2,2)) cy=(akwy+cc*amat(1,2))/(cc*amat(2,2)) return end ; pro cff1n, Bl, bxp, byp, position = position ; s=size(Bl) nx = s(1) ny = s(2) nx2 = 2*nx ny2 = 2*ny alx = nx2/(2.0*!pi) aly = ny2/(2.0*!pi) bl0 = total(Bl)/(nx*ny) rtemp = fltarr(nx2,ny2) + bl0 rtemp( 0:(nx-1), 0:(ny-1) ) = Bl cftbz = fft( rtemp, -1) cx = cftbz cy = cx amat=fltarr(3,3) if keyword_set(position) then begin amat = matrix_tn( position ) endif else begin amat(*,0)= [ 1.0, 0.0, 0.0 ] amat(*,1)= [ 0.0, 1.0, 0.0 ] amat(*,2)= [ 0.0, 0.0, 1.0 ] endelse bn0 = bl0/amat(2,2) bx0 = amat(0,2)*bn0 by0 = amat(1,2)*bn0 for j=1, ny2 do begin for i=1, nx2 do begin if i eq 1 and j eq 1 then goto, endloop kwx = i-1 kwy = j-1 if kwx gt nx then kwx = kwx -nx2 if kwy gt ny then kwy = kwy -ny2 setbtp, kwx,kwy, alx,aly, cx1, cy1, amat cx(i-1,j-1) = cx1 cy(i-1,j-1) = cy1 endloop: endfor endfor kwy = ny for i=1, nx2 do begin kwx = i-1 if kwx gt nx then kwx = kwx -nx2 setbtp, kwx, kwy, alx,aly, cx1, cy1, amat setbtp, kwx,-kwy, alx,aly, cx2, cy2, amat cx(i-1, kwy) = 0.5*(cx1+cx2) cy(i-1, kwy) = 0.5*(cy1+cy2) endfor kwx = nx for j=1, ny2 do begin kwy = j-1 if kwy gt ny then kwy = kwy -ny2 setbtp, kwx,kwy, alx,aly, cx1, cy1, amat setbtp, -kwx,kwy, alx,aly, cx2, cy2, amat cx(kwx, j-1) = 0.5*(cx1+cx2) cy(kwx, j-1) = 0.5*(cy1+cy2) endfor kwx = nx kwy = ny setbtp, kwx, kwy, alx,aly, cx1, cy1, amat setbtp, -kwx, kwy, alx,aly, cx2, cy2, amat setbtp, kwx,-kwy, alx,aly, cx3, cy3, amat setbtp, -kwx,-kwy, alx,aly, cx4, cy4, amat cx(kwx, kwy) = 0.25*(cx1+cx2+cx3+cx4) cy(kwx, kwy) = 0.25*(cy1+cy2+cy3+cy4) cx(0,0) = complex( 0.0, 0.0 ) cy(0,0) = complex( 0.0, 0.0 ) cx = cx*cftbz cy = cy*cftbz cftbz = fft( cx, 1) bxp = float( cftbz( 0:(nx-1), 0:(ny-1) ) ) +bx0 cftbz = fft( cy, 1) byp = float( cftbz( 0:(nx-1), 0:(ny-1) ) ) +by0 print, '(proc.cff1n) end' return end