PRO GET_BLOS, imcube, spcube, line_center=line_center, cont_pos=cont_pos, delta=delta, glande=glande, verbose=verbose, spectfile=spectfile, outdir=outdir, blos_get=blos_get, bmag_get=bmag_get IF (N_PARAMS() LT 2) THEN BEGIN PRINT,'Syntax: GET_BLOS, imcube, spcube, LINE_CENTER=line_center, '+$ 'CONT_POS=cont_pos, DELTA=delta, GLANDE=glande, VERBOSE=verbose, '+$ 'SPECTFILE=spectfile, OUTDIR=outdir, BLOS_GET=blos_get, BMAG_GET=bmag_get' RETURN ENDIF ns = 4 IF (N_ELEMENTS(NLP) NE 1) THEN nlp = 1 IF (N_ELEMENTS(DELTA) NE 1) THEN delta = 1 IF (N_ELEMENTS(CONT_POS) NE 1) THEN cont_pos = 0 IF (N_ELEMENTS(OUTDIR) NE 1) THEN outdir = './' IF (N_ELEMENTS(BLOS_GET) NE 1) THEN blos_get = 1 ;RR lp_readheader,imcube,nx=nx,ny=ny,nt=imnt IF (N_ELEMENTS(SPECTFILE)) THEN BEGIN RESTORE,spectfile lc = WHERE(norm_spect[*,0] EQ MIN(norm_spect[*,0])) nlp = N_ELEMENTS(spect_pos) lp = spect_pos ENDIF ELSE BEGIN lc = line_center[0] wavc = line_center[1] dwav = line_center[2] nlp = line_center[3] lp = (FINDGEN(nlp)-lc) * dwav + wavc ; fact = lp - wavc ; v_dop = c_speed * (lp/lp[lc]-1) ENDELSE l_0 = lp[lc] lp_sel = WHERE(lp NE lp[cont_pos]) lp = lp[lp_sel] nt = imnt / nlp / ns IF KEYWORD_SET(VERBOSE) THEN WINDOW,0 IF KEYWORD_SET(BLOS_GET) THEN blos_cube = FLTARR(nx,ny,nt) IF KEYWORD_SET(BMAG_GET) THEN bmag_cube = FLTARR(nx,ny,nt) im_i_cont = FLTARR(nx,ny,nt) FOR t=0,nt-1 DO im_i_cont[0,0,t] = LP_GET(imcube,t*nlp*ns + cont_pos) pass = 0L totpasses = LONG(nx)*LONG(ny)*LONG(nt) t0 = SYSTIME(/SECONDS) FOR y=0,ny-1 DO BEGIN FOR x=0,nx-1 DO BEGIN spcubepos = LONG(y)*LONG(nx)*LONG(ns) + LONG(x)*LONG(ns) sp_i = LP_GET(spcube,spcubepos) sp_q = LP_GET(spcube,spcubepos + 1L) sp_u = LP_GET(spcube,spcubepos + 2L) sp_v = LP_GET(spcube,spcubepos + 3L) sp_i = FLOAT(sp_i[lp_sel,*]) sp_q = FLOAT(sp_q[lp_sel,*]) sp_u = FLOAT(sp_u[lp_sel,*]) sp_v = FLOAT(sp_v[lp_sel,*]) FOR t=0,nt-1 DO BEGIN ; Determine inclination gamma sp_v_min = MIN(sp_v[*,t],where_sp_v_min, MAX=sp_v_max, SUBSCRIPT_MAX=where_sp_v_max) ; Get minimum and maximum of Stokes V profile ; where_sp_v_min = WHERE(sp_v[*,t] EQ sp_v_min) ; where_sp_v_max = WHERE(sp_v[*,t] EQ sp_v_max) lowbound_p = (where_sp_v_max - delta) > 0 ; Determine boundaries uppbound_p = (where_sp_v_max + delta) < (nlp-2) lowbound_m = (where_sp_v_min - delta) > 0 uppbound_m = (where_sp_v_min + delta) < (nlp-2) lowbound_p = lowbound_p[0] uppbound_p = uppbound_p[0] lowbound_m = lowbound_m[0] uppbound_m = uppbound_m[0] IF KEYWORD_SET(BMAG_GET) THEN BEGIN gamma_vals = FLTARR(4*delta+2) FOR ll=lowbound_p,uppbound_p DO BEGIN xvar_p = sp_v[ll,t] / (SQRT(sp_q[ll,t]^2 + sp_u[ll,t]^2)) ; IF (xvar_p NE 0.) THEN gamma_vals[ll-lowbound_p] = ACOS((SQRT(1+xvar_p^2) - 1) / FLOAT(xvar_p)) ELSE gamma_vals[ll-lowbound_p] = 0. gamma_vals[ll-lowbound_p] = ACOS((SQRT(1+xvar_p^2) - 1) / FLOAT(xvar_p)) ENDFOR FOR ll=lowbound_m,uppbound_m DO BEGIN xvar_m = sp_v[ll,t] / (SQRT(sp_q[ll,t]^2 + sp_u[ll,t]^2)) ; IF (xvar_m NE 0.) THEN gamma_vals[ll-lowbound_m+2*delta+1] = ACOS((SQRT(1+xvar_m^2) - 1) / FLOAT(xvar_m)) ELSE gamma_vals[ll-lowbound_m+2*delta+1] = 0. gamma_vals[ll-lowbound_m+2*delta+1] = ACOS((SQRT(1+xvar_m^2) - 1) / FLOAT(xvar_m)) ENDFOR ; meangamma = MEAN(gamma_vals) meangamma = MEAN(gamma_vals[WHERE(FINITE(gamma_vals) EQ 1)]) IF (FINITE(meangamma) EQ 0) THEN STOP ENDIF ELSE meangamma = 'N/A' ; Determine centroids of right and left circularly polarized components denominator_p = REPLICATE(im_i_cont[x,y,t],nlp-1) - (sp_i[*,t] + sp_v[*,t]) ; stop numerator_p = lp * denominator_p denom_integral_p = TRAPEZ(lp,denominator_p) numer_integral_p = TRAPEZ(lp,numerator_p) l_p = numer_integral_p/denom_integral_p ; stop denominator_m = REPLICATE(im_i_cont[x,y,t],nlp-1) - (sp_i[*,t] - sp_v[*,t]) numerator_m = lp * denominator_m denom_integral_m = TRAPEZ(lp,denominator_m) numer_integral_m = TRAPEZ(lp,numerator_m) l_m = numer_integral_m/denom_integral_m ; stop IF ((FINITE(l_p) EQ 0) OR (FINITE(l_m) EQ 0)) THEN STOP blos_val = ((l_p - l_m)/2.)/(4.67D-13 * glande * l_0^2) blos_cube[x,y,t] = blos_val IF KEYWORD_SET(BMAG_GET) THEN BEGIN bmag_val = blos_val / ABS(COS(meangamma)) bmag_cube[x,y,t] = bmag_val ENDIF ELSE bmag_val = 'N/A' ; stop pass += 1L IF KEYWORD_SET(VERBOSE) THEN $ PROCESS_TIMER,pass,totpasses,t0,EXTRA_OUTPUT='(x,y,t)=('+STRTRIM(x,2)+','+STRTRIM(y,2)+','+STRTRIM(t,2)+'). (l_p,l_m,gamma)=('+STRTRIM(l_p,2)+','+STRTRIM(l_m,2)+','+STRTRIM(meangamma,2)+')'+$ ' (B_los,B_mag)='+STRTRIM(blos_val,2)+','+STRTRIM(bmag_val,2)+')' ENDFOR ENDFOR ENDFOR PRINT,'About to write cube(s)... Press .c to continue.' STOP IF KEYWORD_SET(BLOS_GET) THEN LP_WRITE, blos_cube, outdir+'blos-'+strjoin((strsplit(systime(),' ',/extr))[[4,1,2]])+'.fcube' IF KEYWORD_SET(BMAG_GET) THEN LP_WRITE, bmag_cube, outdir+'bmag-'+strjoin((strsplit(systime(),' ',/extr))[[4,1,2]])+'.fcube' END