; NAME :	nonfff
; PURPOSE :
;	generate a simulated magnetogram of a non-force-free field
; CATEGORY :
; 	solar magnetic field computation package (MAGPACK2)
; CALLING SEQUENCE :
;	.run nonfff
; INPUTS :
;	none ( parameters are given in the program )
; OPTIONAL INPUT PARAMETERS : 
;	none
; KEYWORD PARAMETERS :
;	none
; OUTPUTS :
;	creates a magnetogram file
; COMMON BLOCKS : 	none
; SIDE EFFECTS : 	none
; RESTRICTIONS :	none
; PROCEDURE :
;
; MODIFICATION HISTORY :
;	T.Sakurai, 1992
;
;------------------------------------------------------------------------
;
;  vector potential A=(0,Ay,0), Ay=a0*exp(-ax*x**2-ay*y**2-akz*z)
;  xx: point where the field bb is evaluated
;  xxs: reference position of the field
;
	pro nonfff,  xx,bb,aj, a0,ax,ay,akz,xxs
;
        x=xx(0)-xxs(0)
        y=xx(1)-xxs(1)
        z=xx(2)
        e=exp(-ax*x*x-ay*y*y-akz*z)
	bb = fltarr(3)
	aj = fltarr(3)

        bb(0)= a0*akz*e
        bb(1)= 0.0
        bb(2)=-2.0*a0*ax*x*e
        aj(0)= 4.0*a0*ax*ay*x*y*e
        aj(1)=(2.0*ax-4.0*(ax*x)^2-akz*akz)*a0*e
        aj(2)= 2.0*a0*ay*akz*y*e
        return
	end
;
;	generate a simulated magnetogram based on a non-force-free field
;	main program
;
	bmax = 1000.0
	xxs = [0.0, 0.0, -1.0]
	wx = 1.0
	wy = 1.0
	hz = 1.0

	nx = 21
	ny = 21
	xs = -5.0
	ys = -5.0
	dx = 0.5
	dy = 0.5

	mparam = defmp()

;	mparam.kindm = 23L
	mparam.kindm = 3L
	mparam.nbytem = 2
	mparam.nframe = 4
	mparam.tsys='GMT+09'
	mparam.inst='DUMMY'
	mparam.mfile='NONFFF'
	mparam.nx = nx
	mparam.ny = ny
	mparam.xs = xs
	mparam.ys = ys
	mparam.dx = dx
	mparam.dy = dy
	mparam.mesh = 0
	mparam.title='NONFFF'
;
	rdata = fltarr(nx,ny,8)
	kz=0
	kx=1
	ky=2
	kj = 3
	mparam.status.bl = kz
	mparam.status.bx = kx
	mparam.status.by = ky
	mparam.status.jn = kj
	mparam.factors.bl=0.2
	mparam.factors.bx=0.2
	mparam.factors.by=0.2
	mparam.factors.jn=1.0
;
        a0=1.0
        ax=1.0/wx^2
        ay=1.0/wy^2
        akz=1.0/hz
        bmax1=0.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 )
;
		nonfff,  xx,bb,aj, a0,ax,ay,akz,xxs
;
		rdata(i,j,kx)=bb(0)
		rdata(i,j,ky)=bb(1)
		rdata(i,j,kz)=bb(2)
        	ajn= total( aj*vn )
		rdata(i,j,kj)=ajn
	endfor
	endfor
	bmax1 = sqrt( max( rdata(*,*,0)^2 + rdata(*,*,1)^2 + rdata(*,*,2)^2 ) )
;
	a0 = bmax*a0/bmax1
	rdata(*,*,0:3) = a0*rdata(*,*,0:3)
	if mparam.kindm eq 3L then mparam.status.jn = -1
;
	writem, mparam, rdata, filename='~/idl/nonfff.dat'
end