pro aia_mag_loopsim,fileset,fov,ds,hmax,nr
;+
; Project     : AIA
;
; Name        : AIA_MAG_SIM
;
; Category    : SDO/HMI simulations 
;
; Explanation : simulates current-carrying loops and loop field lines
;
; Syntax      : IDL> aia_mag_loopsim,fileset,fov,ds,hmax
;
; Inputs      : fileset	  = filename
;		fov	  = field-of-view (in solar radii)
;		ds	  =
;               hmax      = 
;
; Outputs     :  
;
; History     : 11-Aug-2011, Version 1 written by Markus J. Aschwanden
;
; Contact     : aschwanden@lmsal.com
;-

;reads loop coefficients________________________________________________
coeff_file=fileset+'_coeff.dat'
readcol,coeff_file,m1,x1,y1,z1
coeff	=[[m1],[x1],[y1],[z1]]

;simulates loop file____________________________________________________
nloop	=nr
npmax	=200
seed1	=1
seed2	=2
ran1	=randomu(seed1,nloop)
ran2	=randomu(seed2,nloop)
loopfile=fileset+'_loop.dat'
simfile =fileset+'_sim.dat'
field_=fltarr(npmax,nloop,4)
openw,2,loopfile
openw,3,simfile

for il=0,nloop-1 do begin
 xf1	=x1[il]
 yf1	=y1[il]
 zf1	=sqrt(1.-xf1^2-yf1^2)
 alpha0 =0.
 aia_mag_fieldline,xf1,yf1,zf1,ds,hmax,coeff,curr0,alpha0,curr,nc   ;-->curr
 xm 	=x1[il]+0.01
 ym 	=y1[il]+0.01
 zm	=sqrt(1.-xm^2-ym^2)
 alpha  =ran2[il]				;units per solar radii
 coeff_powell=[alpha,xf1,yf1]
 print,'POWELL=',coeff_powell 
 printf,3,coeff_powell
 aia_mag_fieldline,xm,ym,zm,ds,hmax,coeff,curr,alpha,field,ns0
 is1	=long(ns0*0.1)
 is2	=long(ns0*(1.+ran1[il])/2.)
 ns	=is2-is1
 ind 	=is1+findgen(ns) 
 for ip=0,ns-1 do begin
  x	=field(is1+ip,0)
  y	=field(is1+ip,1)
  z	=field(is1+ip,2)
  b	=field(is1+ip,3)
  printf,2,il,x,y,z,b
  field_(ip,il,0)=x
  field_(ip,il,1)=y
  field_(ip,il,2)=z
  field_(ip,il,3)=b
 endfor
endfor
close,2
close,3
print,'file written : ',loopfile
print,'Simulation parameters written : ',simfile

;control display________________________________________________________
x1	=fov(0)
x2	=fov(2)
y1	=fov(1)
y2	=fov(3)
z1	=min(z)
z2	=max(z)
grid	=5
rsun	=1.
window,0,xsize=512,ysize=1024
loadct,5
col	=50+200*(1+findgen(nloop))/float(nloop+1)
!x.range=[x1,x2]
!y.range=[y1,y2]
!x.style=1
!y.style=1
!p.position=[0.05,0.05,0.95,0.45]
plot,curr(*,0),curr(*,1),thick=5
for iloop=0,nloop-1 do begin
 ind	=where(field_(*,iloop,2) ne 0) 
 oplot,field_(ind,iloop,0),field_(ind,iloop,1),color=col(iloop),thick=3
endfor
coordsphere2,rsun,0,0,0,0,grid,0,0
!noeras=1

!p.position=[0.05,0.55,0.95,0.95]
!y.range=[z1,z2]+0.5*(x2-x1)*[-1,1]
plot,curr(*,0),curr(*,2),thick=5
for iloop=0,nloop-1 do begin
 ind	=where(field_(*,iloop,2) ne 0) 
 oplot,field_(ind,iloop,0),field_(ind,iloop,2),color=col(iloop),thick=3
endfor
coordsphere2,rsun,0,0,-90,0,grid,0,0
end