code: def compute_operation(file1, file2, operation, output, display,
                        args_z1z2, args_bbox, args_keystitle, args_geometry):
      with fits.open(file1) as hdulist:
          image_header1 = hdulist[0].header
          image1 = hdulist[0].data.astype(np.float)
      naxis1 = image_header1['naxis1']
      naxis2 = image_header1['naxis2']
      if display == 'all':
          ximshow_file(file1.name,
                       args_z1z2=args_z1z2, args_bbox=args_bbox,
                       args_keystitle=args_keystitle,
                       args_geometry=args_geometry,
                       debugplot=12)
      with fits.open(file2) as hdulist:
          image_header2 = hdulist[0].header
          image2 = hdulist[0].data.astype(np.float)
      naxis1_ = image_header2['naxis1']
      naxis2_ = image_header2['naxis2']
      if display == 'all':
          ximshow_file(file2.name,
                       args_z1z2=args_z1z2, args_bbox=args_bbox,
                       args_keystitle=args_keystitle,
                       args_geometry=args_geometry,
                       debugplot=12)
      if naxis1 != naxis1_:
          raise ValueError("NAXIS1 values are different.")
      if naxis2 != naxis2_:
          raise ValueError("NAXIS2 values are different.")
      if operation == "+":
          solution = image1 + image2
      elif operation == "-":
          solution = image1 - image2
      elif operation == "*":
          solution = image1 * image2
      elif operation == "/":
          solution = image1 / image2
      else:
          raise ValueError("Unexpected operation=" + str(operation))
      hdu = fits.PrimaryHDU(solution.astype(np.float), image_header1)
      hdu.writeto(output, overwrite=True)
      if display in ['all', 'result']:
          ximshow_file(output.name,
                       args_z1z2=args_z1z2, args_bbox=args_bbox,
                       args_keystitle=args_keystitle,
                       args_geometry=args_geometry,
                       debugplot=12)