;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Generate the initial constant alpha fields. This begins ; ; by cutting out a section of an HMI FITS file. We then run ; ; a batch file. it is easier to get the aia images with the ; ; mg, just to make sure everything is in place ; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;get mg_unrot by using the program hmi_archive.pro in dana/mpole/anal ;;;;;date is a string with the date of the form 'yyyy-mm-ddThh:mm:ss' mg_unrot=mdi_archives(date) ;;;;;use read_aia() to decompose the .fits files downloaded from the internet AIA171 = read_aia('HMI_fits_file.fits',/mg,/zoom) ;;;;;it has you select an area, it is fine if the aia image is bigger than the mg_unrot image ;;;;; we only use it find loops anyways ;;;;;coalign the images using 1600A, ;;;;;plot the 1600A and then use contour like this to allign contour,mg_unrot.img,mg_unrot.x,mg_unrot.y,level=[-75.0,75.0],/overplot ;;;;;keep track of the dx and dy that you make bacuase you need to change the aia ;;;;;wavelengths that you got, and also store the xshift and yshift in the aia's themselves ;;;;;for some reason I have needed to redefine the aia images once I got the x and yshifts ;;;;;I did it using this code himg=aia335 aia335= {img:himg.img,x:himg.x,y:himg.y,rad:himg.rad,b0:himg.b0,$ p:0.0000,date:himg.date,$ wavelen:himg.wavelen,source:himg.source,exptime:himg.corrected_for_exptime,$ xshift:dx,yshift:dy} ;!;!;!;!; don't forget to suply your own xshift and yshift ;;;;;save all your aia and mg_unrot data to a sav file, so you can restore it later ;;;;;after you have a magnetogram that is unrotated and saved then put it through a batch file like at ;;;;; /nfs/hl2/data/jjensen/ar11112/make_ca_fields.batch ;;;;; I just ran the batch file normal in the IDL window ;;;;; also notice the congrid program, the bigger the number that you divide by, ;;;;; the faster the program runs, but it also pixelates it a lot, so I have kept it to 2-3, ;;;;; and it usually runs in about 6-8 hours, which is just about right, ;;;;; if you make it too long it will affect all the processes afterwards ;;;;; and make the whole process take two to three times as long ;;;;; but it will be a 'more accurate' calculation, but not worth the effort, I think ;;;;; if you make it too short it will not be as accurate a calculation ;; Get 16 million colors instead of 256, makes everything look perty device, decomposed = 1 clr = rainbow() ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Now, we want to trace loops in our AIA images. This is most easily ; ; accomplished by passing the images though a highpass filter. The ; ; highpass.pro function is located in /home/tarrl/work/v_current/. ; ; The function xrt_mdi_job_umka.pro, located in ; ; /home/tarrl/work/lfff/full/ is where we actually trace the loops ; ; using a Bezier spline. Trace until you can't....hence the 'i+=1' line ; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;I found it useful to create a script file that I could run every time I ;;;;;wanted to do these, an example of a the scripts that I wrote for each wavelength ;;;;;can be found at ;;;;; /nfs/hl2/data/jjensen/ar11112/2012_07_30_folder_171_script.txt ;;;;; it contains these steps, and makes it so you can start drawing lines without having ;;;;; to put in all this stuff every time ;; make an aia image run through a highpass filter with something like: himg = aia171[0] himg = {img:float(himg.img),x:himg.x,y:himg.y,rad:himg.rad,b0:himg.b0,$ p:himg.p,date:himg.date,time:himg.time,$ wavelen:himg.wavelen,source:himg.source,exptime:himg.exptime,$ xshift:himg.xshift,yshift:himg.yshift} img = highpass(himg.img, sigma = 5, width = 10) ; highpass lives in /home/tarrl/work/gen_scripts himg.img=img ;; hand find loops with bezier splines: i = -1 xshift = himg.xshift yshift = himg.yshift ;;;;;repeat this next line as many times as you can, ;;;;;and keep an eye on your i or you might overwrite some of your lines i+=1 & xrt_mdi_job_umka, himg, inmodel, outfolder=outfolder, xshift=xshift, yshift=yshift, i, /showprev ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; After you have as many loops as you want, fit each loop to the ; ; various lfffs generated above. I found that the time that was ; ;required depended on your alpha fields that were generated, with ; ;a run of 6-8 hours on the alpha fields thiss process would only ; ;take about 7-10 minutes per line. ; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; files = file_search(outfolder+'fl_???.sav', count=nf) for j = 0, nf-1 do begin $ lothian_n_browning_visual_fit_v2, outfolder, j, d_ave, /nodisp,$ mag_dir='/nfs/hl2/data/jjensen/ar11112/ca_fields/' ;;;;;notice the specification of the mag_dir, this needs to point to the ;;;;; current magnetogram, it is a new thing that was put in the programs ;;;;; in my /home/jjensen/IDL folder, so if you are running the programs ;;;;; from there you need to specify the mag_dir, it pops up in other ;;;;;places too ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Select the best \alpha-h value for each loop, BEING SURE TO USE ; ; THE CORRECT AIA IMAGE AND XSHIFT/YSHIFT VALUES. Selection ; ; algorithms are provided in Anny Malanushenko's dissertation. ; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; for j = 0,nf-1 do begin $ select_minima, aia335[0], outfolder, j, xshift=xshift, yshift=yshift, $ mag_dir='/nfs/hl2/data/jjensen/ar11112/ca_fields/' ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Now, calculate energies. Each step is described in detail in the following. ; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; After selecting the best match for all the loops, ;; calculate the integrals and along each loop dana_energy_paper_prep, '', outfolder, '/nfs/hl2/data/jjensen/ar11112/all_loops/loops.txt','', $ '/nfs/hl2/data/jjensen/ar11112/all_loops/energy.sav',$ mag_dir='/nfs/hl2/data/jjensen/ar11112/ca_fields/' ;; ~/.../energy.sav has a bunch of stuff in it, but mostly we want loops_bdl and loops_bpdl, ;; the integrals along the loops. Our energy estimate is given by ;; \sum_i (\Phi_i \Delta C_i )/(8\pi) ;; where \Delta C_i is (loops_bdl[i] - loops_bpdl[i]) and \Phi_i is found ;; using the voronoi tesselations of the footpoints. ;; First change the magnetograms to centimeters k=6.96e10/mg.rad mg.x=mg.x*k & mg.y=mg.y*k & mg.z=mg.z*k & mg0=mg ;;restore proper energy file that was generated by the dana_enerfy_paper_prep program restore,'/nfs/hl2/data/jjensen/ar11112/all_loops/energy.sav' ;; next, get the voronoi fluxes phi_n = voronoi_phi( loops_x[0,*], loops_y[0,*], mg, /neg, exf=2, rad=0.5,/graph ) phi_p = voronoi_phi( loops_x[1,*], loops_y[1,*], mg, /pos , exf=2, rad=0.5,/graph ) ; now sum! erg_free = total((loops_bdl-loops_bpdl)*(phi_p-phi_n)/2/(8*3.14159D0)) erg_tot = total((loops_bdl)*(phi_p-phi_n)/2/(8*3.14159D0)) N=float(n_elements(loops_bdl)) erg_free_stuff = total(((loops_bdl-loops_bpdl)*(phi_p-phi_n)/2/(8*3.14159D0))^2) erg_tot_stuff = total(((loops_bdl)*(phi_p-phi_n)/2/(8*3.14159D0))^2) erg_free_unc=((N/(N-1))*erg_free_stuff-(erg_free^2/(N-1)))^0.5 erg_tot_unc=((N/(N-1))*erg_tot_stuff-(erg_tot^2/(N-1)))^0.5 help,erg_free,erg_free_unc,erg_tot,erg_tot_unc,n ;; Can also display loops and footpoints: display, mg.img, mg.x, mg.y,/aspect, title = mg.date ;; footpoints: for j = 0,150 do begin $ plots, loops_x[0,j], loops_y[0,j],psym = 1, color =clr[5] & $ plots, loops_x[1,j], loops_y[1,j],psym = 2, color =clr[0] & $ endfor ;; loops: for j = 0,150 do begin $ ll = loops_all[*, loops_ij[0, j]:loops_ij[1, j]] & $ plots, ll[0,*], ll[1,*], ll[2,*], color = clr[5] & $ endfor ;;;;;;;;;;;;;;;; viewing transforms for over plotting stuff ;;;;;;;;;;;;;;; vnew = solar_rotate_view(varr[113], to=aia.date) fl_view_xform, cline1,vnew fl_view_xform, cline3,vnew show_fp, spls[113],nls113 ,/nosepx,/over,view=vnew fl_view_xform, cline2,vne ;;;;;this is a graph that I generated with the four separate wavelengths ;;;;; I will leave it here for reference ;;;;;;;;;;;;;;;;;;;;;;;;;;ploting the final graph of energy;;;;;;;;;;;;;;;;;;;;; ;;;;;initialize the graph with the proper data points ;set up for graphing nice data this is for ar11166 no open loops erg_free=[2.42e+32,3.66e+31,3.18e+32,6.88e+31] erg_free_unc=[1.59e+32,1.86e+31,1.69e+32,3.26e+31] erg_tot=[9.72e+32,7.15e+32,9.77e+32,6.36e+32] erg_tot_unc=[2.73e+32,1.80e+32,2.58e+32,1.42e+32] all_free=[2.51e+32,2.51e+32,2.51e+32,2.51e+32] all_free_unc=[1.06e+32,1.06e+32,1.06e+32,1.06e+32] all_tot=[9.62e+32,9.62e+32,9.62e+32,9.62e+32] all_tot_unc=[1.73e+32,1.73e+32,1.73e+32,1.73e+32] ;initialize data for ar11092 no open loops 0500 UT erg_free=[1.41e+32,8.38e+31,1.15e+32,1.925e+32] erg_free_unc=[1.05e+32,3.18e+31,4.10e+31,1.14e+32] erg_tot=[4.21e+32,5.14e+32,4.73e+32,4.98e+32] erg_tot_unc=[1.40e+32,9.40e+31,8.88e+31,1.60e+32] all_free=[2.59e+32,2.59e+32,2.59e+32,2.59e+32] all_free_unc=[9.84e+31,9.84e+31,9.84e+31,9.84e+31] all_tot=[7.21e+32,7.21e+32,7.21e+32,7.21e+32] all_tot_unc=[1.40e+32,1.40e+32,1.40e+32,1.40e+32] ;initialize data for ar110922 at 1700 UT july 1 erg_free=[5.58e+31,1.73e+31,7.32e+31,1.17e+32] erg_free_unc=[3.19e+31,6.42e+30,2.84e+31,6.01e+31] erg_tot=[3.66e+32,2.70e+32,3.89e+32,3.97e+32] erg_tot_unc=[5.82e+31,4.18e+31,5.60e+31,1.13e+32] all_free=[1.00e+32,1.00e+32,1.00e+32,1.00e+32] all_free_unc=[2.84e+31,2.84e+31,2.84e+31,2.84e+31] all_tot=[4.72e+32,4.72e+32,4.72e+32,4.72e+32] all_tot_unc=[5.05e+31,5.05e+31,5.05e+31,5.05e+31] ;initialize data for ar11112 16 Oct 2010 16:50 UT erg_free=[4.71e+30,8.44e+30,4.63e+30,1.72e+31] erg_free_unc=[2.73e+30,3.67e+30,1.74e+30,1.26e+31] erg_tot=[6.56e+31,7.10e+31,6.29e+31,6.29e+31] erg_tot_unc=[1.55e+31,1.48e+31,1.27e+31,2.61e+31] all_free=[1.13e+31,1.13e+31,1.13e+31,1.13e+31] all_free_unc=[4.68e+30,4.68e+30,4.68e+30,4.68e+30] all_tot=[8.45e+31,8.45e+31,8.45e+31,8.45e+31] all_tot_unc=[1.23e+31,1.23e+31,1.23e+31,1.23e+31] ;this is the FREE energy plot x=[0,1,2,3] plot,erg_free,charsize=1.75,psym=5,yrange=[1.0e+30,4.0e+32],xrange=[-.5,3.5],$ xtitle='171 angstrom 193 angstrom 211 angstrom 335 angstrom',$ ytitle='ergs',title='Calculated Free Energies for ar11092 0500 UT Aug 1',xminor=1 oploterr, x, erg_free, +erg_free_unc, psym=3,/HIBAR,thick=2 oploterr, x, erg_free, +erg_free_unc, psym=3,/LOBAR,thick=2 oplot,all_free,linestyle=5,thick=3 oplot,all_free-all_free_unc,linestyle=1,thick=3 oplot,all_free+all_free_unc,linestyle=1,thick=3 ;this is the TOTAL energy plot x=[0,1,2,3] plot,erg_tot,charsize=1.75,psym=5,yrange=[3.0e+31,1.2e+32],xrange=[-.5,3.5],$ xtitle='171 angstrom 193 angstrom 211 angstrom 335 angstrom',$ ytitle='ergs',title='Calculated Total Energies for ar11112 1650 UT oct 16',xminor=1 oploterr, x, erg_tot, +erg_tot_unc, psym=3,/HIBAR,thick=2 oploterr, x, erg_tot, +erg_tot_unc, psym=3,/LOBAR,thick=2 oplot,all_tot,linestyle=5,thick=3 oplot,all_tot-all_tot_unc,linestyle=1,thick=3 oplot,all_tot+all_tot_unc,linestyle=1,thick=3 ;;;;;This is a little snipit of code that I used to step through the lines that I ;;;;;had drawn and see if they were up to par ;;;;;;;;;;;;;;;;;;;;fixing the open loop stuff;;;;;;;;;;;;;;;;;;;;;;;;;;; i+=1 restore,'fl_'+string(i, format='(I3.3)')+'_proc_man.sav' dx=mg.x[1]-mg.x[0] & dy=mg.y[1]-mg.y[0] & dz=mg.z[1]-mg.z[0] b_=best_line_man & b_[0, *]=b_[0, *]*dx+mg.x[0] & $ b_[1, *]=b_[1, *]*dy+mg.y[0] & b_[2, *]=b_[2, *]*dz+mg.z[0] display,mg.img,mg.x,mg.y,/aspect oplot, b_[0, *], b_[1, *], color=255, lines=2, thick=2 print,i