################################################################################
# Combine a coronal solution with a heliospheric solution.
################################################################################

# The directory that contains the (r,t,p) meshes to use for
# the combined solution.
set mesh_dir=MeshCoarse

set dir_corona=../corona_corotating_frame
set dir_helio=../heliosphere_shifted_3au

set seq_corona=002
set seq_helio=002

# Tag for the output file names.
set root="C${seq_corona}+H${seq_helio}"

# Radius at which to join the two solutions.
set r_join=27.5

# List of fields to process.
set fields=(br bt bp vr vt vp p rho t)

################################################################################
# Begin processing.
################################################################################

set mesh="-x ${mesh_dir}/r.hdf -y ${mesh_dir}/t.hdf -z ${mesh_dir}/p.hdf"

# Make a dummy data set with the correct structure.

remesh -oval 0 ${mesh} ${dir_corona}/br${seq_corona}.hdf tmp_r.hdf

# Set up masks for the coronal and helio solutions.

sds_func -func one -x tmp_r.hdf tmp_r.hdf
sds_func -func scale_x -x tmp_r.hdf tmp_r.hdf

sds_func -func mask_le -a ${r_join} -x tmp_r.hdf tmp_MC.hdf
sds_func -func ax+b -a -1 -b 1 -x tmp_MC.hdf tmp_MH.hdf

\rm -f tmp_r.hdf

# Interpolate the two solutions.

foreach fld (${fields})
  remesh -oval 0 ${mesh} ${dir_corona}/${fld}${seq_corona}.hdf tmp_${fld}1.hdf
  remesh -oval 0 ${mesh} ${dir_helio}/${fld}${seq_helio}.hdf tmp_${fld}2.hdf
  sds_func -func 'ax*y+b' -a 1 -b 0 -x tmp_${fld}1.hdf -y tmp_MC.hdf tmp_${fld}1.hdf
  sds_func -func 'ax*y+b' -a 1 -b 0 -x tmp_${fld}2.hdf -y tmp_MH.hdf tmp_${fld}2.hdf
  sds_func -func ax+by -a 1 -b 1 -x tmp_${fld}1.hdf -y tmp_${fld}2.hdf ${fld}${root}.hdf
  \rm -f tmp_${fld}1.hdf tmp_${fld}2.hdf
end

\rm -f tmp_*.hdf
