;====================================
pro save_img, data, sx,sy,folder,file
;====================================
siz=size(data)
print,'Writing ',folder,file
print,'Columns= ',siz[1],' Rows= ',siz[2]
print,'Min= ',min(data),' Max= ',max(data)
device,decomposed=0
loadct,39
tvlct,rr,gg,bb,/get
window,xsize=sx,ysize=sy,/free
tvscl,data
img24=tvrd(true=1)
img8=color_quan(img24,1,rr,gg,bb)
write_png,folder+file+'.png',img8,rr,gg,bb
openw, lunw, folder+file+'.txt',/get_lun,width=5000
printf,lunw, data
close,lunw
free_lun, lunw
openw, lunw, folder+file+'.bin',/get_lun
writeu,lunw, data
close,lunw
free_lun, lunw
end

;=================
function vec_mag,v
;=================
return, sqrt(v[0]^2+v[1]^2+v[2]^2)
end

;=================
pro magnetic_field
;=================
n=100.0       ;Number of loops in each coil
current=2.0 ; current in amps
radius=0.20 ; Radius of coil
a=0.10      ;Tube to center of coil distance
cons=1.26e-6/(4.0*!PI)
dtheta=2.0*!PI/1000.0
xylim=100
bx=replicate(0.0,2*xylim+1,2*xylim+1)
by=replicate(0.0,2*xylim+1,2*xylim+1)
bz=replicate(0.0,2*xylim+1,2*xylim+1)
for xp=-xylim, xylim do begin
  print,'xp= ',xp
  for yp=-xylim, xylim do begin
    rp=float([xp,yp,0.0])/1000.0
    for icoil=0,1 do begin ; Do both coils
      if (icoil eq 0) then aa=-a else aa=a
      ra=[0.0,0.0,aa]
      for itheta=0,1000 do begin ; Integrate around each coil (2PI)
        theta=(float(itheta)/1000.0)*2.0*!PI
        ri=[radius*cos(theta),radius*sin(theta),0.0]+ra
        dl=[-radius*dtheta*sin(theta),radius*dtheta*cos(theta),0.0]
        r=rp-ri
        rmag=vec_mag(r)
        ;Biot and Savart Law
        db=(cons*current*n*crossp(dl,r)/rmag^3)
        bx[xp+xylim,yp+xylim]=bx[xp+xylim,yp+xylim]+db[0]
        by[xp+xylim,yp+xylim]=by[xp+xylim,yp+xylim]+db[1]
        bz[xp+xylim,yp+xylim]=bz[xp+xylim,yp+xylim]+db[2]
      endfor
    endfor
  endfor
endfor


folder='c:\Documents and Settings\henry.snyder\physics\labs\magfield\'
file='bx'
save_img, bx,2*xylim+1,2*xylim+1,folder,file
file='by'
save_img, by,2*xylim+1,2*xylim+1,folder,file
file='bz'
save_img, bz,2*xylim+1,2*xylim+1,folder,file
end