Add return_footprint choice to function reproject_and_coadd

Version: 0.14.0
```
1 coaddfile = os.path.join(path, 'new_cube.fits')
2 header_out = hdus[0].header
3 array = reproject_and_coadd(hdus,
4 wcs_out,
5 #output_array=hdu_out[0].data.astype(np.float64),
6 shape_out=shape_out,
7 combine_function = 'mean',
8 reproject_function=reproject_interp,
9 match_background=True,
10 return_footprint=False,
11 progress_bar = tqdm,
12 block_size=(2048, 2048),
13 parallel=cpu_count(),
14 roundtrip_coords=False,)
15 fits.writeto(coaddfile, array.astype(np.float32), header=hdu_out[0].header, overwrite=True)

File ~/.local/lib/python3.10/site-packages/reproject/mosaicking/coadd.py:314, in reproject_and_coadd(input_data, output_projection, shape_out, input_weights, hdu_in, hdu_weights, reproject_function, combine_function, match_background, background_reference, output_array, output_footprint, block_sizes, progress_bar, blank_pixel_value, intermediate_memmap, **kwargs)
310 array = footprint = None
312 logger.info(f"Calling {reproject_function.__name__} with shape_out={shape_out_indiv}")
314 array, footprint = reproject_function(
315 (array_in, wcs_in),
316 output_projection=wcs_out_indiv,
317 shape_out=shape_out_indiv,
318 hdu_in=hdu_in,
319 output_array=array,
320 output_footprint=footprint,
321 **kwargs,
322 )
324 if weights_in is not None:
326 if intermediate_memmap:

ValueError: too many values to unpack (expected 2)
```

Start in reproject/mosaicking/coadd.py around line 314, where reproject_function returns array and footprint, and reproduce the reported call from the issue with return_footprint=False. Trace how reproject_and_coadd handles that return value; done means the requested return_footprint choice no longer causes the reported unpacking error and produces the intended return shape.

Written by the indexing model from the issue text.