Estimating the flare location coordinates

[paolomassa]

I think it could be nice to have a little wrapper around the code for estimating the flare location, which is subsequently used to compute the visibilities. Below a simple example, which could be the starting point for implementing an appropriate function

def estimate_flare_location(cpd, time_range, 
                            isc = [3, 20, 22, 16, 14, 32, 21, 26, 4, 24, 8, 28],
                            energy_range = [6,10],
                            imsize = [512, 512]*u.pixel,
                            make_plot = True):
    
    # Create calibrated visibilities
    meta_pixels = create_meta_pixels(cpd, time_range=time_range,
                                     energy_range=energy_range, phase_center=[0, 0] * u.arcsec, 
                                     no_shadowing=True)
    
    vis = create_visibility(meta_pixels)
    
    uu, vv = get_visibility_info_giordano()
    vis.u = uu
    vis.v = vv
    cal_vis = calibrate_visibility(vis)
    
    idx = np.argwhere(np.isin(cal_vis.isc, isc)).ravel()

    stix_vis = Visibility(vis=cal_vis.obsvis[idx], u=cal_vis.u[idx], v=cal_vis.v[idx])
    skeys = stix_vis.__dict__.keys()
    for k, v in cal_vis.__dict__.items():
        if k not in skeys:
            setattr(stix_vis, k, v[idx])
    
    # Compute Back Projection
    rsun_obs = cpd.meta['rsun_arc']*u.arcsec
    pixel = rsun_obs * 2.6 / imsize.value
    
    fd_bp_map = vis_to_map(stix_vis, imsize, pixel_size=pixel)
    fd_bp_map.meta['rsun_obs'] = rsun_obs.value
    
    max_pixel = np.argwhere(fd_bp_map.data == fd_bp_map.data.max()).ravel() *u.pixel
    
    # because WCS axes are reverse order
    stx_flare_loc = fd_bp_map.pixel_to_world(max_pixel[1], max_pixel[0])
    
    if make_plot:
        fig = plt.figure()
        axes = fig.add_subplot(projection=fd_bp_map)
        fd_bp_map.plot(axes=axes)
        fd_bp_map.draw_limb()
        axes.plot_coord(stx_flare_loc, marker='.', markersize=50, fillstyle='none', color='r', markeredgewidth=2)
        axes.set_xlabel('STIX X')
        axes.set_ylabel('STIX Y')
    
    return stx_flare_loc`
    
Maybe part of this function could be implemented in xrayvision as general function for determining flare location from visibilities