code: import os
  import sys
  import numpy as np
  import matplotlib.pyplot as plt
  from astropy.coordinates import SkyCoord
  from astropy import units as u
  
  sys.path.append("/usr/custom/pyLIMA-1.0.0")
  from pyLIMA import event, telescopes, microlmodels
  import MulensModel as mm
  
  
  # Common settings:
  t_0 = 2456900.
  u_0 = 0.1
  t_E = 150.
  pi_E_N = 1.
  pi_E_E = 2.
  ra_deg = 270.
  dec_deg = -30.
  time = np.linspace(t_0-80., t_0+80., 1000)
  d_time = 2450000.
  
  # MulensModel calculations and plot:
  coords = SkyCoord(ra_deg, dec_deg, unit=u.deg)
  params = {'t_0': t_0, 'u_0': u_0, 't_E': t_E,
            'pi_E_N': pi_E_N, 'pi_E_E': pi_E_E}
  mm_model = mm.Model(params, coords=coords)
  mm_mag = mm_model.magnification(time)
  plt.plot(time-d_time, mm_mag, 'r.', label='MulensModel')
  
  # Read VBBL output and plot it:
  vbbl_data = np.loadtxt("fake.out", unpack=True)
  plt.plot(vbbl_data[0], vbbl_data[1], 'g-.', label='VBBL')
  
  # This are the changes I have to make to make the results as close as possible:
  pi_E_E = -pi_E_E
  pi_E_N = -pi_E_N
  
  # pyLIMA calculations and plots:
  your_event = event.Event()
  your_event.ra = ra_deg
  your_event.dec = dec_deg
  data_1 = np.array([time, time*0.+15., time*0.+.01]).T
  telescope_1 = telescopes.Telescope(light_curve_magnitude=data_1)
  your_event.telescopes.append(telescope_1)
  your_event.check_event()
  model_1 = microlmodels.create_model(
      'PSPL', your_event, parallax=['Annual', t_0])
  model_1.define_model_parameters()
  pyLIMA_parameters = model_1.compute_pyLIMA_parameters(
      [t_0, u_0, t_E, pi_E_N, pi_E_E])
  model = model_1.compute_the_microlensing_model(telescope_1, pyLIMA_parameters)
  mag_pyLIMA = model_1.model_magnification(telescope_1, pyLIMA_parameters)
  plt.plot(time-d_time, mag_pyLIMA, 'b.', label='pyLIMA')
  
  # Compare pyLIMA and MM:
  index = np.argmax(np.abs(mm_mag - mag_pyLIMA))
  print("Largest difference is for: ", index, time[index]-d_time)
  print("pyLIMA:", mag_pyLIMA[index])
  print("MM:", mm_mag[index])
  
  # This is the end:
  plt.legend(loc='best')
  plt.xlabel('JD-2450000')
  plt.ylabel('magnification')
  plt.show()
   """
  Calculates interpolation tables for elliptical integral of
  the first and second kind.
  """
  import math
  import numpy as np
  from math import sin, cos, sqrt
  from scipy import integrate
  from scipy.interpolate import interp1d
  from scipy.special import ellipk, ellipe
  # These are complete elliptic integrals of the first and the second kind.
  
  
  accuracy = 1.e-6
  n_divide = 10 + 1
  x_start = 2.e-5
  x_stop = 1. - 1.e-12
  n_start = 10
  
  # Settings end here.
  
  
  def get_ellip(x):
      k = []
      e = []
      for x_ in x:
          if x_ not in get_ellip.k:
              get_ellip.k[x_] = ellipk(x_)
              get_ellip.e[x_] = ellipe(x_)
          k.append(get_ellip.k[x_])
          e.append(get_ellip.e[x_])
      return (np.array(k), np.array(e))
  get_ellip.k = dict()
  get_ellip.e = dict()
  
  x = np.logspace(np.log10(x_start), np.log10(x_stop), n_start)
  
  iteration = 0
  add = [None]
  while len(add) > 0:
      iteration += 1
      add = []
      (k, e) = get_ellip(x)
      interp_k = interp1d(np.log10(x), k, kind='cubic')
      interp_e = interp1d(np.log10(x), e, kind='cubic')
      for i in range(len(x)-1):
          x_1 = x[i]
          x_2 = x[i+1]
          check = np.logspace(np.log10(x_1), np.log10(x_2), n_divide)[1:-1]
          (check_true_k, check_true_e) = get_ellip(check)
          check_k = []
          check_e = []
          for c in check:
              check_k.append(interp_k(np.log10(c)))
              check_e.append(interp_e(np.log10(c)))
          check_k = np.array(check_k)
          check_e = np.array(check_e)
          relative_diff_k = np.abs(check_k - check_true_k) / check_true_k
          relative_diff_e = np.abs(check_e - check_true_e) / check_true_e
          relative_diff_max = np.amax(
              np.array([relative_diff_k, relative_diff_e]), axis=0)
          index = np.argsort(relative_diff_max)[-1]
          if relative_diff_max[index] < accuracy:
              continue
          add.append(check[index])
      new_x = np.sort(add + x.tolist())
      x = new_x
  
  for (x_, k_, e_) in zip(x, k, e):
      print(x_, k_, e_)
  
apis: matplotlib.pyplot.legend numpy.abs matplotlib.pyplot.show matplotlib.pyplot.ylabel numpy.array matplotlib.pyplot.plot numpy.linspace matplotlib.pyplot.xlabel numpy.loadtxt
apis: scipy.special.ellipk scipy.special.ellipe numpy.abs numpy.argsort numpy.log10 numpy.array