Skip to content
  • Search
Skins
  • Light
  • Cerulean
  • Cosmo
  • Flatly
  • Journal
  • Litera
  • Lumen
  • Lux
  • Materia
  • Minty
  • Morph
  • Pulse
  • Sandstone
  • Simplex
  • Sketchy
  • Spacelab
  • United
  • Yeti
  • Zephyr
  • Dark
  • Cyborg
  • Darkly
  • Quartz
  • Slate
  • Solar
  • Superhero
  • Vapor

  • Default (No Skin)
  • No Skin
Collapse

ZMT zurich med tech

  1. Home
  2. Sim4Life
  3. Analysis & Postprocessing
  4. Max modulation tool logic

Max modulation tool logic

Scheduled Pinned Locked Moved Analysis & Postprocessing
4 Posts 3 Posters 723 Views
  • Oldest to Newest
  • Newest to Oldest
  • Most Votes
Reply
  • Reply as topic
Log in to reply
This topic has been deleted. Only users with topic management privileges can see it.
  • D Offline
    D Offline
    dstoupis
    wrote on last edited by dstoupis
    #1

    How does the Max Modulation tool calculate the maximum modulation of two electric fields? I am asking since it is not described in the manual and it also does not appear in any of the tutorials.
    Based on the this post, the tool is based on the formula derived in the Grossman et al., 2017 paper which is the following:

    render.png

    I applied this formula and the results I got were not correct when post-processing the two electric fields in Python, outside of Sim4Life. The above formula is for the |E_1| > |E_2| area. The code that I am using in Python is:

    import numpy as np
    import scipy.io as io
    
    base_data = io.loadmat('Directory of extracted mat file for the simulation 1')
    delta_data = io.loadmat('Directory of extracted mat file for the simulation 2')
    
    e_field_1 = np.nan_to_num(base_data['EM_E_x_y_z_f0_Snapshot0'])
    e_field_2 = np.nan_to_num(delta_data['EM_E_x_y_z_f0_Snapshot0'])
    
    # Calculate the subtraction and addtion of the vectors
    E_plus = np.add(e_field_1, e_field_2) # Create the sum of the E fields
    E_minus = np.subtract(e_field_1, e_field_2) # Create the difference of the E fields
    
    # Calculate the angles between the two fields for each vector
    # Assuming that E_a is the difference of the vectors
    dot_angle = np.einsum('ij,ij->i', e_field_1, E_minus)
    cross_angle = np.linalg.norm(np.cross(e_field_1, E_minus), axis=1)
    
    angles = np.arctan2(cross_angle, dot_angle)
    
    # Condition to have two times the E2 field amplitude
    max_condition = np.linalg.norm(e_field_2, axis=1) < np.linalg.norm(e_field_1, axis=1)*np.cos(angles)
    
    e1_gr_e2 = np.where(np.linalg.norm(e_field_1, axis=1) > np.linalg.norm(e_field_2, axis=1), max_condition, False)
    
    # Calculate the maximum envelope
    envelope = np.where(e1_gr_e2, 2.0*np.linalg.norm(e_field_2, axis=1), envelope)
    
    # Find the other points where the condition is not true
    tmp = np.logical_xor(e1_gr_e2, np.linalg.norm(e_field_1, axis=1) > np.linalg.norm(e_field_2, axis=1))
    envelope = np.where(tmp, 2.0*(np.linalg.norm(np.cross(e_field_2, E_minus), axis=1)/np.linalg.norm(E_minus, axis=1)), envelope)
    
    np.nan_to_num(envelope)
    

    The resulting graphs can be seen in the figure below, with the calculations from the above code shown on the left-most image, results from Sim4Life shown in the middle and the common points between Sim4Life and Python calculations shown in the right-most image.

    SimAndClalc2.png

    I have tried many variations of the code, setting the amplitudes of electric field on and two to double the value for all points and trying to find the common points and also tried different variations with the indices in the cross product formula. Still with all the variations I could not get all the points filled in the common points plot. The image with all calculated common points is shown below. The correct simulation image is also shown below.

    Sec_req.png

    Thank you!

    1 Reply Last reply
    1
    • D Offline
      D Offline
      dstoupis
      wrote on last edited by dstoupis
      #2

      Well.... after long time finally the solution has been found! Based on the Rampersad et al. 2019 where the mathematics behind the tTIS is better explained, the key line is that any of the two fields has to change direction whenever the angle between the two vectors is not less than 90 degrees. That was the magic line! Below I post the fully functional code as a python function. The arguments e_field_1 and e_field_2 are n, 3 arrays of the electric fields as produced by Sim4Life.

      def modulation_envelope(e_field_1, e_field_2, dir_vector=None):
          E_plus = np.add(e_field_1, e_field_2) # Create the sum of the E fields
          
          if dir_vector is None:
              envelope = np.zeros(e_field_1.shape[0])
      
              # Calculate the angles between the two fields for each vector
              dot_angle = np.einsum('ij,ij->i', e_field_1, e_field_2)
              cross_angle = np.linalg.norm(np.cross(e_field_1, e_field_2), axis=1)
              angles = np.arctan2(cross_angle, dot_angle)
              
              # Flip the direction of the electric field if the angle between the two is greater or equal to 90 degrees
              e_field_2 = np.where(np.broadcast_to(angles >= np.pi/2., (3, e_field_2.shape[0])).T, -e_field_2, e_field_2)
              
              # Recalculate the angles
              dot_angle = np.einsum('ij,ij->i', e_field_1, e_field_2)
              cross_angle = np.linalg.norm(np.cross(e_field_1, e_field_2), axis=1)
              angles = np.arctan2(cross_angle, dot_angle)
      
              E_minus = np.subtract(e_field_1, e_field_2) # Create the difference of the E fields
              
              # Condition to have two times the E2 field amplitude
              max_condition = np.linalg.norm(e_field_2, axis=1) < np.linalg.norm(e_field_1, axis=1)*angs
              e1_gr_e2 = np.where(np.linalg.norm(e_field_1, axis=1) > np.linalg.norm(e_field_2, axis=1), max_condition, False)
      
              # Condition to have two times the E1 field amplitude
              max_condition_2 = np.linalg.norm(e_field_1, axis=1) < np.linalg.norm(e_field_2, axis=1)*angs
              e2_gr_e1 = np.where(np.linalg.norm(e_field_2, axis=1) > np.linalg.norm(e_field_1, axis=1), max_condition_2, False)
              
              # Double magnitudes
              envelope = np.where(e1_gr_e2, 2.0*np.linalg.norm(e_field_2, axis=1), envelope) # 2E2 (First case)
              envelope = np.where(e2_gr_e1, 2.0*np.linalg.norm(e_field_1, axis=1), envelope) # 2E1 (Second case)
              
              # Cross product
              # Redifine e1_gr_e2 to match the area ouside where the max_condition is met
              e1_gr_e2 = np.where(np.linalg.norm(e_field_1, axis=1) > np.linalg.norm(e_field_2, axis=1), np.logical_not(max_condition), False)
              envelope = np.where(e1_gr_e2, 2.0*(np.linalg.norm(np.cross(e_field_2, E_minus), axis=1)/np.linalg.norm(E_minus, axis=1)), envelope) # (First case)
              
              e2_gr_e1 = np.where(np.linalg.norm(e_field_2, axis=1) > np.linalg.norm(e_field_1, axis=1), np.logical_not(max_condition_2), False)
              envelope = np.where(e2_gr_e1, 2.0*(np.linalg.norm(np.cross(-e_field_1, -E_minus), axis=1)/np.linalg.norm(-E_minus, axis=1)), envelope) # (Second case)
          else:
              envelope = np.abs(np.abs(np.dot(E_plus, dir_vector)) - np.abs(np.dot(E_minus, dir_vector)))
          return envelope
      
      1 Reply Last reply
      2
      • SylvainS Offline
        SylvainS Offline
        Sylvain
        ZMT
        wrote on last edited by
        #3

        Thanks a lot for sharing your solution!

        1 Reply Last reply
        0
        • L Offline
          L Offline
          LJ
          wrote on last edited by
          #4
          This post is deleted!
          1 Reply Last reply
          0
          Reply
          • Reply as topic
          Log in to reply
          • Oldest to Newest
          • Newest to Oldest
          • Most Votes


          • Login

          • Don't have an account? Register

          • Login or register to search.
          • First post
            Last post
          0
          • Search