; NAME :	dipole
; PURPOSE :
;	generate a simulated magnetogram of a dipole field
; CATEGORY :
; 	solar magnetic field computation package (MAGPACK2)
; CALLING SEQUENCE :
;	.run dipole
; INPUTS :
;	none (parameters are specified in the program)
; OPTIONAL INPUT PARAMETERS : 
;	none
; KEYWORD PARAMETERS :
;	none
; OUTPUTS :
;	data are written to a file
; COMMON BLOCKS : 	none
; SIDE EFFECTS : 	none
; RESTRICTIONS :	none
; PROCEDURE :
;
; MODIFICATION HISTORY :
;	T.Sakurai, 1992
;
;------------------------------------------------------------------------
;
;  dipole field
;  xx: point where the field bb is evaluated
;  xxs: source point
;  dmom: dipole moment
;  ddir: direction of dipole moment (unit vector)
;
	pro dipole,  xx,bb,dmom,xxs,ddir
;
	dp1 = fltarr(3)
	dp2 = fltarr(3)
;
	if abs(ddir(0)) lt abs(ddir(1))  and $ 
		abs(ddir(0)) lt abs(ddir(2))  then begin
		dp1(1)=1.0
	endif else begin
	if abs(ddir(1)) lt abs(ddir(0))  and $ 
		abs(ddir(1)) lt abs(ddir(2))  then begin
		dp1(2)=1.0
	endif else begin
		dp1(2)=1.0
	endelse
	endelse
;
	dp2(0)=ddir(1)*dp1(2)-ddir(2)*dp1(1)
	dp2(1)=ddir(2)*dp1(0)-ddir(0)*dp1(2)
	dp2(2)=ddir(0)*dp1(1)-ddir(1)*dp1(0)
;
	dp2 = dp2/ sqrt( total(dp2*dp2) )
;
	dp1(0)=dp2(1)*ddir(2)-dp2(2)*ddir(1)
	dp1(1)=dp2(2)*ddir(0)-dp2(0)*ddir(2)
	dp1(2)=dp2(0)*ddir(1)-dp2(1)*ddir(0)
;
	xx0 = xx - xxs
	r=sqrt( total( xx0*xx0 ) )
	x=xx0(0)*dp1 (0)+xx0(1)*dp1 (1)+xx0(2)*dp1 (2)
	y=xx0(0)*dp2 (0)+xx0(1)*dp2 (1)+xx0(2)*dp2 (2)
	z=xx0(0)*ddir(0)+xx0(1)*ddir(1)+xx0(2)*ddir(2)
	bx=dmom*3.0*x*z/r^5
	by=dmom*3.0*y*z/r^5
	bz=dmom*(2.0*z*z-x*x-y*y)/r^5
;
	bb = bx*dp1 + by*dp2 + bz*ddir
	return
	end
;
;	generate a simulated magnetogram based on a dipole field
;	main program
;
	bmax = 1000.0
	ddir = [1.0, 0.0, 0.0]
	xxs = [0.0, 0.0, -1.0]
	nx = 21
	ny = 21
	xs = -5.0
	ys = -5.0
	dx = 0.5
	dy = 0.5

	mparam = defmp()

;	mparam.kindm = 2L
	mparam.kindm = 3L
	mparam.nbytem = 2
	mparam.nframe = 3
	mparam.tsys='GMT+09'
	mparam.inst='DUMMY'
	mparam.mfile='DIPOLE'
	mparam.nx = nx
	mparam.ny = ny
	mparam.xs = xs
	mparam.ys = ys
	mparam.dx = dx
	mparam.dy = dy
	mparam.mesh = 0
	mparam.title='DIPOLE'
;
	rdata = fltarr(nx,ny,8)
	kz=0
	kx=1
	ky=2
	mparam.status.bl = kz
	mparam.status.bx = kx
	mparam.status.by = ky
	mparam.factors.bl=0.2
	mparam.factors.bx=0.2
	mparam.factors.by=0.2
;
	ddir = ddir/ sqrt( total(ddir*ddir) )
;
	dmom=1.0

	for j=0,ny-1 do begin
	for i=0,nx-1 do begin
		xx = [ xs+dx*i, ys+dy*j, 0.0 ]
		xx = projt( xx, mparam, vn )
;
		dipole, xx,bb,dmom,xxs,ddir
;
		rdata(i,j,kx)=bb(0)
		rdata(i,j,ky)=bb(1)
		rdata(i,j,kz)=bb(2)
	endfor
	endfor
	bmax1 = sqrt( max( rdata(*,*,0)^2 + rdata(*,*,1)^2 + rdata(*,*,2)^2 ) )
;
	dmom=bmax/bmax1
	rdata(*,*,0:2) = dmom*rdata(*,*,0:2)
;
	writem, mparam, rdata, filename='~/idl/dipole.dat'
end