Skip to content

Commit

Permalink
Fix issue in VISFWDFIT_PSO update for upper and lower bound (#229)
Browse files Browse the repository at this point in the history 
* Update for upper and lower bound
* Update comments
  • Loading branch information
@annavolp
annavolp authored Nov 28, 2024
1 parent 035af0f commit c5a3310
Show file tree
Hide file tree
Showing 7 changed files with 185 additions and 77 deletions.
95 changes: 78 additions & 17 deletions stix/idl/processing/imaging/stx_vis_fwdfit_pso.pro
View file
Original file line number Diff line number Diff line change
Expand Up @@ -78,6 +78,8 @@ function stx_vis_fwdfit_pso, configuration, vis, aux_data, $

n_sources = n_elements(configuration)

flare_loc_HPC = stx_hpc2stx_coord(vis[0].xyoffset, aux_data, /inverse) ;Helioprojective Cartesian coordinate

if keyword_set(SRCIN) then begin

if n_circle gt 0 then begin
Expand Down Expand Up @@ -108,23 +110,21 @@ function stx_vis_fwdfit_pso, configuration, vis, aux_data, $
; Set up file I/O error handling.
ON_IOError, y_stix_c
; Cause type conversion error.
if flag2 then this_y_c = double(srcin.circle[j].param_opt.param_y)

if flag2 then this_y_c = double(srcin.circle[j].param_opt.param_y)

if (((flag1 eq 0) and (flag2 eq 1)) or ((flag1 eq 1 ) and (flag2 eq 0))) then begin
Catch, /Cancel
message, "Fix both x and y positions or none of them."
endif
if ((flag1 eq 1) and (flag2 eq 1)) then begin

;From HPC to STIX coordinate frame
this_xy_c = stx_hpc2stx_coord([double(srcin.circle[j].param_opt.param_x),double(srcin.circle[j].param_opt.param_y)], aux_data)

srcin.circle[j].param_opt.param_x = string(this_xy_c[0])
srcin.circle[j].param_opt.param_y = string(this_xy_c[1])

endif


endif

endfor
endif

Expand Down Expand Up @@ -162,14 +162,13 @@ function stx_vis_fwdfit_pso, configuration, vis, aux_data, $
message, "Fix both x and y positions or none of them."
endif
if ((flag3 eq 1) and (flag4 eq 1)) then begin

;From HPC to STIX coordinate frame
this_xy_e = stx_hpc2stx_coord([double(srcin.ellipse[j].param_opt.param_x),double(srcin.ellipse[j].param_opt.param_y)], aux_data)

srcin.ellipse[j].param_opt.param_x = string(this_xy_e[0])
srcin.ellipse[j].param_opt.param_y = string(this_xy_e[1])

endif


flag=1
Catch, theError
Expand Down Expand Up @@ -221,15 +220,13 @@ function stx_vis_fwdfit_pso, configuration, vis, aux_data, $
message, "Fix both x and y positions or none of them."
endif
if ((flag5 eq 1) and (flag6 eq 1)) then begin

;From HPC to STIX coordinate frame
this_xy_l = stx_hpc2stx_coord([double(srcin.loop[j].param_opt.param_x),double(srcin.loop[j].param_opt.param_y)], aux_data)

srcin.loop[j].param_opt.param_x = string(this_xy_l[0])
srcin.loop[j].param_opt.param_y = string(this_xy_l[1])

endif


endif

flag=1
Catch, theError
Expand All @@ -249,12 +246,76 @@ function stx_vis_fwdfit_pso, configuration, vis, aux_data, $

endif else begin

srcin = vis_fwdfit_pso_multiple_src_create(vis, configuration)
srcin = vis_fwdfit_pso_multiple_src_create(vis, configuration, flare_loc_HPC)

endelse

if n_circle gt 0 then begin

for j=0, n_circle-1 do begin
;the solution for the position is sought in a rectangle centered in [0.,0.]
srcin.circle[j].lower_bound.l_b_x = srcin.circle[j].lower_bound.l_b_x - flare_loc_HPC[0]
srcin.circle[j].upper_bound.u_b_x = srcin.circle[j].upper_bound.u_b_x - flare_loc_HPC[0]
srcin.circle[j].lower_bound.l_b_y = srcin.circle[j].lower_bound.l_b_y - flare_loc_HPC[1]
srcin.circle[j].upper_bound.u_b_y = srcin.circle[j].upper_bound.u_b_y - flare_loc_HPC[1]

; upper and lower bound transformation: from the Solar Orbiter coordinate frame to the STIX coordinate frame
this_l_b_x = srcin.circle[j].lower_bound.l_b_y
this_u_b_x = srcin.circle[j].upper_bound.u_b_y
this_l_b_y = - srcin.circle[j].upper_bound.u_b_x
this_u_b_y = - srcin.circle[j].lower_bound.l_b_x

srcin.circle[j].lower_bound.l_b_x = this_l_b_x
srcin.circle[j].upper_bound.u_b_x = this_u_b_x
srcin.circle[j].lower_bound.l_b_y = this_l_b_y
srcin.circle[j].upper_bound.u_b_y = this_u_b_y
endfor
endif

if n_ellipse gt 0 then begin
for j=0, n_ellipse-1 do begin
;the solution for the position is sought in a rectangle centered in [0.,0.]
srcin.ellipse[j].lower_bound.l_b_x = srcin.ellipse[j].lower_bound.l_b_x - flare_loc_HPC[0]
srcin.ellipse[j].upper_bound.u_b_x = srcin.ellipse[j].upper_bound.u_b_x - flare_loc_HPC[0]
srcin.ellipse[j].lower_bound.l_b_y = srcin.ellipse[j].lower_bound.l_b_y - flare_loc_HPC[1]
srcin.ellipse[j].upper_bound.u_b_y = srcin.ellipse[j].upper_bound.u_b_y - flare_loc_HPC[1]

; upper and lower bound transformation: from the Solar Orbiter coordinate frame to the STIX coordinate frame
this_l_b_x = srcin.ellipse[j].lower_bound.l_b_y
this_u_b_x = srcin.ellipse[j].upper_bound.u_b_y
this_l_b_y = - srcin.ellipse[j].upper_bound.u_b_x
this_u_b_y = - srcin.ellipse[j].lower_bound.l_b_x

srcin.ellipse[j].lower_bound.l_b_x = this_l_b_x
srcin.ellipse[j].upper_bound.u_b_x = this_u_b_x
srcin.ellipse[j].lower_bound.l_b_y = this_l_b_y
srcin.ellipse[j].upper_bound.u_b_y = this_u_b_y
endfor
endif


if n_loop gt 0 then begin
for j=0, n_loop-1 do begin
;the solution for the position is sought in a rectangle centered in [0.,0.]
srcin.loop[j].lower_bound.l_b_x = srcin.loop[j].lower_bound.l_b_x - flare_loc_HPC[0]
srcin.loop[j].upper_bound.u_b_x = srcin.loop[j].upper_bound.u_b_x - flare_loc_HPC[0]
srcin.loop[j].lower_bound.l_b_y = srcin.loop[j].lower_bound.l_b_y - flare_loc_HPC[1]
srcin.loop[j].upper_bound.u_b_y = srcin.loop[j].upper_bound.u_b_y - flare_loc_HPC[1]

; upper and lower bound transformation: from the Solar Orbiter coordinate frame to the STIX coordinate frame
this_l_b_x = srcin.loop[j].lower_bound.l_b_y
this_u_b_x = srcin.loop[j].upper_bound.u_b_y
this_l_b_y = - srcin.loop[j].upper_bound.u_b_x
this_u_b_y = - srcin.loop[j].lower_bound.l_b_x

srcin.loop[j].lower_bound.l_b_x = this_l_b_x
srcin.loop[j].upper_bound.u_b_x = this_u_b_x
srcin.loop[j].lower_bound.l_b_y = this_l_b_y
srcin.loop[j].upper_bound.u_b_y = this_u_b_y
endfor
endif

param_out = vis_fwdfit_pso(configuration, this_vis, srcin, $
param_out = vis_fwdfit_pso(configuration, this_vis, srcin, aux_data.ROLL_ANGLE * !dtor, flare_loc_HPC, $
n_birds = n_birds, tolerance = tolerance, maxiter = maxiter, $
uncertainty = uncertainty, $
imsize=imsize, pixel=pixel, $
Expand All @@ -278,7 +339,7 @@ function stx_vis_fwdfit_pso, configuration, vis, aux_data, $

nsrc = N_ELEMENTS(srcstr)
FOR n = 0, nsrc-1 DO BEGIN

; From STIX coordinate frame to HPC one
xy_hpc = stx_hpc2stx_coord([srcstr[n].srcx,srcstr[n].srcy], aux_data, /inverse)

srcstr[n].srcx = xy_hpc[0]
Expand All @@ -298,7 +359,7 @@ function stx_vis_fwdfit_pso, configuration, vis, aux_data, $
srcstr[n].srcpa, $
srcstr[n].srcx, srcstr[n].srcy, srcstr[n].loop_angle]
PRINT, n+1, srcstr[n].srctype, temp, FORMAT="(I5, A13, F13.2, 1F13.1, F12.1, 2F11.1, F11.1, 2F12.1)"

;; Note that the position uncertainty is exchanged due to the 90 degree rotation. (see stx_solo2stx_coord.pro for more details)
temp = [ fitsigmas[n].srcflux,fitsigmas[n].srcfwhm_max, fitsigmas[n].srcfwhm_min, $
fitsigmas[n].srcpa, fitsigmas[n].srcy, fitsigmas[n].srcx, fitsigmas[n].loop_angle]
PRINT, ' ', '(std)', temp, FORMAT="(A7, A11, F13.2, 1F13.1, F12.1, 2F11.1, F11.1, 2F12.1)"
Expand Down
49 changes: 30 additions & 19 deletions stix/idl/processing/imaging/vis_fwdfit_func_pso.pro
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,7 @@ function vis_fwdfit_func_pso, xx, extra = extra
configuration = extra.configuration
param_opt = extra.param_opt
mapcenter = extra.mapcenter
roll_angle = extra.roll_angle

loc_circle = where(configuration eq 'circle', n_circle)>0
loc_ellipse = where(configuration eq 'ellipse', n_ellipse)>0
Expand Down Expand Up @@ -68,7 +69,9 @@ function vis_fwdfit_func_pso, xx, extra = extra
; Set up file I/O error handling.
ON_IOError, bad_x_circle
; Cause type conversion error.
if flag then xx[*, 4*i+1] = xx[*, 4*i+1] * 0. + double(param_opt[4*i+1]) - mapcenter[0]
; If the x and y coordinate are fixed, take into account the telescope rotation and transform the coordinates into the telescope coordinate frame.
; The co-ordinates must be in a rectangle centered in [0.,0.], so the mapcenter must be subtracted
if flag then xx[*, 4*i+1] = xx[*, 4*i+1] * 0. + cos(roll_angle)*(double(param_opt[4*i+1])-mapcenter[0]) - sin(roll_angle)*(double(param_opt[4*i+2])-mapcenter[1])
x_loc = reform(xx[*, 4*i+1], [n_particles,1])

flag=1
Expand All @@ -81,7 +84,9 @@ function vis_fwdfit_func_pso, xx, extra = extra
; Set up file I/O error handling.
ON_IOError, bad_y_circle
; Cause type conversion error.
if flag then xx[*, 4*i+2] = xx[*, 4*i+2] * 0. + double(param_opt[4*i+2]) - mapcenter[1]
; If the x and y coordinate are fixed, take into account the telescope rotation and transform the coordinates into the telescope coordinate frame.
; The co-ordinates must be in a rectangle centered in [0.,0.], so the mapcenter must be subtracted
if flag then xx[*, 4*i+2] = xx[*, 4*i+2] * 0. + sin(roll_angle)*(double(param_opt[4*i+1])-mapcenter[0]) + cos(roll_angle)*(double(param_opt[4*i+2])-mapcenter[1])
y_loc = reform(xx[*,4*i+2], [n_particles,1])

flag=1
Expand All @@ -97,9 +102,10 @@ function vis_fwdfit_func_pso, xx, extra = extra
if flag then xx[*, 4*i+3] = xx[*, 4*i+3] * 0. + double(param_opt[4*i+3])
fwhm = reform(xx[*,4*i+3], [n_particles,1])

vispred_re += flux * exp(-(!pi^2. * fwhm^2. / (4.*alog(2.)))#(u^2. + v^2.))*cos(2*!pi*((x_loc#u)+(y_loc#v)))
vispred_im += flux * exp(-(!pi^2. * fwhm^2. / (4.*alog(2.)))#(u^2. + v^2.))*sin(2*!pi*((x_loc#u)+(y_loc#v)))

; Take into account the telescope rotation and transform the coordinates into the telescope coordinate frame.
vispred_re += flux * exp(-(!pi^2. * fwhm^2. / (4.*alog(2.)))#(u^2. + v^2.))*cos(2*!pi*(((cos(roll_angle) * x_loc + sin(roll_angle) * y_loc)#u)+((-sin(roll_angle) * x_loc + cos(roll_angle) * y_loc))#v))
vispred_im += flux * exp(-(!pi^2. * fwhm^2. / (4.*alog(2.)))#(u^2. + v^2.))*sin(2*!pi*(((cos(roll_angle) * x_loc + sin(roll_angle) * y_loc)#u)+((-sin(roll_angle) * x_loc + cos(roll_angle) * y_loc))#v))

endfor

endif
Expand Down Expand Up @@ -202,7 +208,9 @@ function vis_fwdfit_func_pso, xx, extra = extra
; Set up file I/O error handling.
ON_IOError, bad_x_ellipse
; Cause type conversion error.
if flag then xx[*, n_circle*4+6*i+1] = xx[*, n_circle*4+6*i+1] * 0. + double(param_opt[n_circle*4+6*i+1]) - mapcenter[0]
; If the x and y coordinate are fixed, take into account the telescope rotation and transform the coordinates into the telescope coordinate frame.
; The co-ordinates must be in a rectangle centered in [0.,0.], so the mapcenter must be subtracted
if flag then xx[*, n_circle*4+6*i+1] = xx[*, n_circle*4+6*i+1] * 0. + cos(roll_angle)*(double(param_opt[n_circle*4+6*i+1])-mapcenter[0]) - sin(roll_angle)*(double(param_opt[n_circle*4+6*i+2])-mapcenter[1])
x_loc = reform(xx[*,n_circle*4+6*i+1], [n_particles,1])

flag=1
Expand All @@ -215,12 +223,15 @@ function vis_fwdfit_func_pso, xx, extra = extra
; Set up file I/O error handling.
ON_IOError, bad_y_ellipse
; Cause type conversion error.
if flag then xx[*, n_circle*4+6*i+2] = xx[*, n_circle*4+6*i+2] * 0. + double(param_opt[n_circle*4+6*i+2]) - mapcenter[1]
; If the x and y coordinate are fixed, take into account the telescope rotation and transform the coordinates into the telescope coordinate frame.
; The co-ordinates must be in a rectangle centered in [0.,0.], so the mapcenter must be subtracted
if flag then xx[*, n_circle*4+6*i+2] = xx[*, n_circle*4+6*i+2] * 0. + sin(roll_angle)*(double(param_opt[n_circle*4+6*i+1])-mapcenter[0]) + cos(roll_angle)*(double(param_opt[n_circle*4+6*i+2])-mapcenter[1])
y_loc = reform(xx[*,n_circle*4+6*i+2], [n_particles,1])

vispred_re += flux * exp(-(!pi^2. / (4.*alog(2.)))*((u1 * fwhmmajor)^2. + (v1 * fwhmminor)^2.))*cos(2*!pi*((x_loc#u)+(y_loc#v)))
vispred_im += flux * exp(-(!pi^2. / (4.*alog(2.)))*((u1 * fwhmmajor)^2. + (v1 * fwhmminor)^2.))*sin(2*!pi*((x_loc#u)+(y_loc#v)))

; Take into account the telescope rotation and transform the coordinates into the telescope coordinate frame.
vispred_re += flux * exp(-(!pi^2. / (4.*alog(2.)))*((u1 * fwhmmajor)^2. + (v1 * fwhmminor)^2.))*cos(2*!pi*((( cos(roll_angle) * x_loc + sin(roll_angle) * y_loc)#u)+( (-sin(roll_angle) * x_loc + cos(roll_angle) * y_loc)#v)))
vispred_im += flux * exp(-(!pi^2. / (4.*alog(2.)))*((u1 * fwhmmajor)^2. + (v1 * fwhmminor)^2.))*sin(2*!pi*(( (cos(roll_angle) * x_loc + sin(roll_angle) * y_loc)#u)+( (-sin(roll_angle) * x_loc + cos(roll_angle) * y_loc)#v)))

endfor

endif
Expand All @@ -242,9 +253,7 @@ function vis_fwdfit_func_pso, xx, extra = extra
; Cause type conversion error.
if flag then xx[*, n_circle*4+n_ellipse*6+7*i] = xx[*, n_circle*4+n_ellipse*6+7*i] * 0. + double(param_opt[n_circle*4+n_ellipse*6+7*i])
flux = xx[*, n_circle*4+n_ellipse*6+7*i]
;flux = reform(flux, [n_particles,1])
;flux = flux # ones


eccos = reform(xx[*,n_circle*4+n_ellipse*6+7*i+4], [n_particles,1])
ecsin = reform(xx[*,n_circle*4+n_ellipse*6+7*i+5], [n_particles,1])

Expand Down Expand Up @@ -295,9 +304,6 @@ function vis_fwdfit_func_pso, xx, extra = extra
Catch, /Cancel
bad_alpha_loop:
pa = atan(ecsin, eccos) * !radeg
; pa = eccen * 0.
; ind = where(eccen GT 0.001)
; pa[ind] = ATAN(ecsin[ind], eccos[ind]) * !RADEG
flag=0
ENDIF
; Set up file I/O error handling.
Expand All @@ -322,7 +328,9 @@ function vis_fwdfit_func_pso, xx, extra = extra
; Set up file I/O error handling.
ON_IOError, bad_x_loop
; Cause type conversion error.
if flag then xx[*, n_circle*4+n_ellipse*6+7*i+1] = xx[*, n_circle*4+n_ellipse*6+7*i+1] * 0. + double(param_opt[n_circle*4+n_ellipse*6+7*i+1]) - mapcenter[0]
; If the x and y coordinate are fixed, take into account the telescope rotation and transform the coordinates into the telescope coordinate frame.
; The co-ordinates must be in a rectangle centered in [0.,0.], so the mapcenter must be subtracted
if flag then xx[*, n_circle*4+n_ellipse*6+7*i+1] = xx[*, n_circle*4+n_ellipse*6+7*i+1] * 0. + cos(roll_angle)*(double(param_opt[n_circle*4+n_ellipse*6+7*i+1])-mapcenter[0]) - sin(roll_angle)*(double(param_opt[n_circle*4+n_ellipse*6+7*i+2])-mapcenter[1])
x_loc = reform(xx[*,n_circle*4+n_ellipse*6+7*i+1], [n_particles,1])

flag=1
Expand All @@ -335,7 +343,9 @@ function vis_fwdfit_func_pso, xx, extra = extra
; Set up file I/O error handling.
ON_IOError, bad_y_loop
; Cause type conversion error.
if flag then xx[*, n_circle*4+n_ellipse*6+7*i+2] = xx[*, n_circle*4+n_ellipse*6+7*i+2] * 0. + double(param_opt[n_circle*4+n_ellipse*6+7*i+2]) - mapcenter[1]
; If the x and y coordinate are fixed, take into account the telescope rotation and transform the coordinates into the telescope coordinate frame.
; The co-ordinates must be in a rectangle centered in [0.,0.], so the mapcenter must be subtracted
if flag then xx[*, n_circle*4+n_ellipse*6+7*i+2] = xx[*, n_circle*4+n_ellipse*6+7*i+2] * 0. + sin(roll_angle)*(double(param_opt[n_circle*4+n_ellipse*6+7*i+1])-mapcenter[0]) + cos(roll_angle)* (double(param_opt[n_circle*4+n_ellipse*6+7*i+2])-mapcenter[1])
y_loc = reform(xx[*,n_circle*4+n_ellipse*6+7*i+2], [n_particles,1])

flag=1
Expand All @@ -350,7 +360,8 @@ function vis_fwdfit_func_pso, xx, extra = extra
; Cause type conversion error.
if flag then xx[*, n_circle*4+n_ellipse*6+7*i+6] = xx[*, n_circle*4+n_ellipse*6+7*i+6] * 0. + double(param_opt[n_circle*4+n_ellipse*6+7*i+6])
loop_angle = reform(xx[*,n_circle*4+n_ellipse*6+7*i+6], [n_particles,1])


; Take into account the telescope rotation and transform the coordinates into the telescope coordinate frame.
vis_pred = vis_fwdfit_pso_func_makealoop( flux, xx[*,n_circle*4+n_ellipse*6+7*i+3], eccen, x_loc, y_loc, pa, loop_angle, u, v)

vispred_re += vis_pred[*,0:n_elements(u)-1]
Expand Down
Loading

0 comments on commit c5a3310

Please sign in to comment.