pro nlfff_fit,dir,savefile,test ;+ ; Project : SOHO/MDI, SDO/HMI, STEREO ; ; Name : NLFFF_FIT ; ; Category : Magnetic data modeling ; ; Explanation : Fitting of twisted field lines (nonlinear force-free ; field) NLFFF-model to observed loop projections. ; ; Syntax : IDL>nlfff_fit,dir,savefile,test ; ; Inputs : dir,savefile,test ; ; Outputs ; updates savefile ; ; History : 16-May-2011, Version 1 written by Markus J. Aschwanden ; : 17-Dec-2015, Successful test free energy step flare 2014-Mar-29 ; ; 2-Jul-2019, new analytical solution NLFFF_VECTOR.PRO ; ; 7-Oct-2019, independent fit in each time interval, no smoothing ; ; 7-Oct-2019, initial condition coeff2(im,4)=0 ; ; 9-Oct-2019, nmag_np_var=iter < nmag_np ; ; 15-Oct-2019, eliminate weight ; ; 8-Nov-2019, add index.nh,nseg,hmin,hmax,nitmin,nitmax ; ; 8-Nov-2019, eliminate grad_norm,s_len,s_fit,xl,yl,zl ; ; 9-Nov-2019, loop geometries with both linear and curved profiles ; ; Contact : aschwanden@lmsal.com ;- t1 =systime(0,/seconds) print,'___________________________NLFFF_FIT_____________________________' restore,dir+savefile ;--> PARA, INPUT instr =input.instr nmag_p =input.nmag_p ;number of magnetic charges nmag_np =input.nmag_np ;number of magnetic charges amis =input.amis nitmin =input.nitmin nitmax =input.nitmax nh =input.nh nseg =input.nseg hmin =input.hmin hmax =input.hmax hrange =(hmax-hmin) ;_________________EXTRACT LOOP COORDINATES AND VECTORS______________ nf =n_elements(ns_loop_det) ;number of loops if (nf eq 0) then goto,save_data ns_loop =ns_loop_det wave_loop =wave_loop_det field_loop=field_loop_det nfit =nf*nseg len =fltarr(nf) ;proj loop length x_fit =fltarr(nfit,3) ;fit spline points v_fit =fltarr(nfit,3) ;fit spline points for k=0,nf-1 do begin ;number of loops nsmax =max(ns_loop) ns =ns_loop(k) dx =field_loop(ns-1,0,k)-field_loop(0,0,k) dy =field_loop(ns-1,0,k)-field_loop(0,0,k) len(k) =sqrt(dx^2+dy^2) for iseg=0,nseg-1 do begin if (nseg ge 2) then im=long(1+(ns-3)*float(iseg)/float(nseg-1)) < (ns-2) if (nseg eq 1) then im=1 ifit =k*nseg+iseg x_fit(ifit,0) =field_loop(im,0,k) x_fit(ifit,1) =field_loop(im,1,k) x_fit(ifit,2) =field_loop(im,2,k) v_fit(ifit,0) =field_loop(im+1,0,k)-field_loop(im-1,0,k) v_fit(ifit,1) =field_loop(im+1,1,k)-field_loop(im-1,1,k) v_fit(ifit,2) =field_loop(im+1,2,k)-field_loop(im-1,2,k) endfor endfor ;________________LOOP GEOMETRIES_____________________________ nc =long(nh) nn =long(nc*nh*(nh-1)) h_curv =fltarr(nn) s_curv =fltarr(nn) r_curv =fltarr(nn) ns =200 ;typical length s_seg =(long(1+(ns-3)*findgen(iseg)/float(nseg-1)) < (ns-2))/float(ns) h_seg_norm=fltarr(nseg,nn) disp =0 eps =1.e-3 ii =0 for ia=0,1 do begin ;conjugate footpoints for ih2=0,nh-1 do begin for ih1=0,ih2-1 do begin qh1 =float(ih1)/float(nh-1) qh2 =float(ih2)/float(nh-1) ;altitude levels for ic=0,nc-1 do begin ;curvature radii th =atan(qh2) ;inclination angle d =sqrt((1./2.)^2+(qh2/2.)^2) r_curv1=-1/sin(th > eps) ;neg curv radius limit r_curv2=+d/cos(th) ;pos curv radius limit q_curv1=1./r_curv1 ;reciprocal value q_curv2=1./r_curv2 ;reciprocal value q_curv =q_curv1+(q_curv2-q_curv1)*float(ic)/float(nc-1) r_curv[ii]=1./q_curv dd =sqrt(r_curv[ii]^2-d^2) sign_curv=r_curv[ii]/abs(r_curv[ii]) s_curv[ii]=1./2. +dd*sin(th)*sign_curv h_curv[ii]=qh2/2.-dd*cos(th)*sign_curv root =r_curv[ii]^2-(s_curv[ii]-s_seg)^2 h_seg0=(h_curv[ii]+sqrt(root)*sign_curv) > 0. h_seg1=qh1+(qh2-qh1)*h_seg0 if (ia eq 0) then h_seg_norm(*,ii)=h_seg1 if (ia eq 1) then h_seg_norm(*,ii)=reverse(h_seg1) ;...........................TEST DISPLAY...................... if (test eq 2) then begin if (disp eq 0) then begin window,0,xsize=512,ysize=512 clearplot loadct,5 plot,[0,1],[0,0],yrange=[0,1],xrange=[0,1] disp=1 endif oplot,s_seg,h_seg_norm(*,ii) endif ;............................................................. ii =ii+1 endfor endfor endfor endfor print,'Number of geometries = ',nn if (test eq 2) then begin read,'Continue? [0,1] : ',yes if (yes eq 0) then STOP endif ;______________________ITERATION START________________________ angle =90. slope =-1. coeff(*,4)=0. ;reset a0 coeff_best=coeff angle_iter=fltarr(nitmax>1) if (nitmax eq 0) then goto,end_iter for iter=0,nitmax-1 do begin ;_______________________ALTITUDE EVALUATION______________________ if (iter eq 0) then begin xx_fit=fltarr(nfit,3,nn) vv_fit=fltarr(nfit,3,nn) for ii=0,nn-1 do begin ;geometries for k=0,nf-1 do begin ;loops ifit =k*nseg+findgen(nseg) i1fit=k*nseg+((findgen(nseg)-1)>0) i2fit=k*nseg+((findgen(nseg)+1)<(nseg-1)) i3 =k*nseg i4 =k*nseg+(nseg-1) h_seg=h_seg_norm(*,ii)*len(k) r_seg=1.0+hmin+h_seg xx_fit[ifit,0,ii]=x_fit[ifit,0] xx_fit[ifit,1,ii]=x_fit[ifit,1] xx_fit[ifit,2,ii]=sqrt(r_seg^2-x_fit[ifit,0]^2-x_fit[ifit,1]^2) vv_fit[ifit,0,ii]=xx_fit[i2fit,0,ii]-xx_fit[i1fit,0,ii] vv_fit[ifit,1,ii]=xx_fit[i2fit,1,ii]-xx_fit[i1fit,1,ii] vv_fit[ifit,2,ii]=xx_fit[i2fit,2,ii]-xx_fit[i1fit,2,ii] vv_fit(i3,2,ii)=2.*(xx_fit(i3+1,2,ii)-xx_fit(i3,2,ii)) vv_fit(i4,2,ii)=2.*(xx_fit(i4,2,ii) -xx_fit(i4-1,2,ii)) endfor endfor endif dev_geo=fltarr(nf)+90. dev_geo_nseg=fltarr(nseg,nf)+90. for ii=0,nn-1 do begin ;looping geometries nlfff_vector_v3,coeff,xx_fit[*,*,ii],bfff_ ;B-field at pos x_fit vector_product_array,vv_fit[*,*,ii],bfff_,v3,v3_norm,a_rad dev_rad=a_rad < (!pi-a_rad) ;180-deg ambiguity dev_deg=(180./!pi)*dev_rad ;radian into degrees for k=0,nf-1 do begin ifit =k*nseg+findgen(nseg) dev_med=median(dev_deg(ifit)) ;median per length if (dev_med lt dev_geo(k)) then begin dev_geo(k)=dev_med dev_geo_nseg(*,k)=dev_deg(ifit) x_fit(ifit,0)=xx_fit[ifit,0,ii] ;update x-coordinate x_fit(ifit,1)=xx_fit[ifit,1,ii] ;update y-coordinate x_fit(ifit,2)=xx_fit[ifit,2,ii] ;update z-coordinate v_fit(ifit,0)=vv_fit[ifit,0,ii] ;update x-vector v_fit(ifit,1)=vv_fit[ifit,1,ii] ;update y-vector v_fit(ifit,2)=vv_fit[ifit,2,ii] ;update z-vector ns =ns_loop[k] xfit=long((nseg-1)*findgen(ns)/float(ns-1)) field_loop(0:ns-1,2,k)=interpol(x_fit(ifit,2),ifit,xfit) endif endfor endfor if (test eq 2) then begin for k=0,nf-1 do begin print,dev_geo_nseg(*,k) print,'median misalignment = ',dev_geo(k) nlfff_fit_disp,x_det,s_det,k,para,coeff endfor endif ;______________MAXIMUM MISALIGNMENT LIMIT_____________________________ if (iter eq 0) then begin ind_good0=where(dev_geo le amis,ngood0) nf0 =nf ;initial loop number endif qiter =(1.-float(iter)/float(nitmin))>0. ngood =long(ngood0+(nf0-ngood0)*qiter) isort =sort(dev_geo) amis_ngood =dev_geo(isort(ngood-1)) ind_good =where(dev_geo le amis_ngood,ngood) ;loop indices nf =ngood ;______________SELF-SELECTION OF GOOD LOOPS____________________________ len_old =len x_old =x_fit v_old =v_fit xx_old =xx_fit vv_old =vv_fit ns_loop =ns_loop(ind_good) wave_loop =wave_loop(ind_good) field_loop=field_loop(*,*,ind_good) len =len_old(ind_good) nfit =nf*nseg x_fit =fltarr(nfit,3) ;fit spline points v_fit =fltarr(nfit,3) ;fit spline points xx_fit =fltarr(nfit,3,nn) ;fit spline points vv_fit =fltarr(nfit,3,nn) ;fit spline points for k=0,nf-1 do begin kold =ind_good(k) iold =nseg*kold+findgen(nseg) ifit =nseg*k +findgen(nseg) x_fit(ifit,*)=x_old(iold,*) v_fit(ifit,*)=v_old(iold,*) xx_fit(ifit,*,*)=xx_old(iold,*,*) vv_fit(ifit,*,*)=vv_old(iold,*,*) endfor ;______________________ALPHA GRADIENT OPTIMIZATION_____________ nmag_np_var=(iter > 2) < nmag_np ;VERSION 5, 2019-Oct-9 grad_alpha=fltarr(nmag_np_var) for im=0,nmag_np_var-1 do begin coeff2=coeff nlfff_vector_v3,coeff2,x_fit,bfff_ ;B-field at position x_fit vector_product_array ,v_fit,bfff_,v3,v3_norm,angle_rad dev_rad =angle_rad < (!pi-angle_rad) ;180-deg ambiguity dev_deg0=(180./!pi)*dev_rad ;radian into degrees dev_med0=median(dev_deg0) grad_alpha[im]=angle-dev_med0 endfor coeff(0:nmag_np_var-1,4)=coeff(0:nmag_np_var-1,4)+grad_alpha(0:nmag_np_var-1) ;______________________NEW_MISALIGNMENT ANGLE_________________ nlfff_vector_v3,coeff,x_fit,bfff_ ;B-field at position x_fit vector_product_array ,v_fit,bfff_,v3,v3_norm,angle_rad dev_rad =angle_rad < (!pi-angle_rad) ;180-deg ambiguity dev_deg1=(180./!pi)*dev_rad ;radian into degrees dev_med1=median(dev_deg1) angle_iter[iter]=dev_med1 if (dev_med1 lt angle) then begin angle =dev_med1 dev_deg_best=dev_deg1 ;array coeff_best=coeff endif ;_______________________CONVERGENCE CRITERION__________________ if (iter ge 1) then begin iter1 =(iter-nitmin-1)>0 iter_ =findgen(nitmax) c =linfit(iter_(iter1:iter),angle_iter(iter1:iter)) slope =c(1) endif print,'ITER, NLOOP, MISALIGN, SLOPE --> ',$ iter,ngood,median(dev_deg_best),slope,format='(A,2I6,2f8.2)' if (slope ge -0.1) and (iter ge nitmin) then goto,end_iter endfor ;for iter=0,nitmax-1 END_ITER: coeff =coeff_best dev_deg =dev_deg_best ;_______________________2D MISALIGNMENT ANGLE__________________ vector_product_array2,v_fit,bfff_,v3,v3_norm,angle_rad2 dev_rad2=angle_rad2 < (!pi-angle_rad2) ;180-deg ambiguity dev_deg2=(180./!pi)*dev_rad2 ;radian into degrees dev_med2=median(dev_deg2) ;median angle angle2 =dev_med2 ;_______________________CPU time_________________________________ t2 =systime(0,/seconds) cpu =t2-t1 print,'FORWARD FITTING : CPU time = ',cpu,format='(a,f8.1)' ;_______________________SAVE PARAMETERS________________________ SAVE_DATA: para.iter =iter para.cpu =cpu para.nitmax =nitmax save,filename=dir+savefile,input,para,bzmap,bzfull,bzmodel,coeff,$ field_loop_det,ns_loop_det,wave_loop_det,$ field_loop,ns_loop,wave_loop,dev_deg,dev_deg2,angle,angle2 print,'parameters saved in file = ',savefile end