From e9987014da2636a3cdb2744ffc1f319acb17c8ce Mon Sep 17 00:00:00 2001
From: cheritier <cheritier@pcp2-cheritier.osupytheas.fr>
Date: Tue, 9 Jan 2024 16:41:26 +0100
Subject: [PATCH] major improvement to use off-axis targets and allow study of
 anisoplanetism

---
 OOPAO/Atmosphere.py | 43 ++++++++++++++++++++++++++++++++-----------
 OOPAO/Telescope.py  | 10 ++++++++--
 2 files changed, 40 insertions(+), 13 deletions(-)

diff --git a/OOPAO/Atmosphere.py b/OOPAO/Atmosphere.py
index 86fb61e..f41165f 100644
--- a/OOPAO/Atmosphere.py
+++ b/OOPAO/Atmosphere.py
@@ -122,6 +122,7 @@ def __init__(self,telescope,r0:float,L0:float,windSpeed:list,fractionalR0:list,w
         self.nExtra                 = 2                 # number of extra pixel to generate the phase screens
         self.wavelength             = 500*1e-9          # Wavelengt used to define the properties of the atmosphere
         self.telescope              = telescope         # associated telescope object
+
         if self.telescope.src is None:
             raise AttributeError('The telescope was not coupled to any source object! Make sure to couple it with an src object using src*tel')
         self.mode                   = mode              # DEBUG -> first phase screen generation mode
@@ -137,6 +138,8 @@ def __init__(self,telescope,r0:float,L0:float,windSpeed:list,fractionalR0:list,w
 # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ATM INITIALIZATION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
     def initializeAtmosphere(self,telescope):
         phase_support = self.initialize_phase_support()
+        self.fov      = telescope.fov
+        self.fov_rad  = telescope.fov_rad
         if self.hasNotBeenInitialized:
             self.initial_r0 = self.r0
             for i_layer in range(self.nLayer):       
@@ -201,10 +204,10 @@ def buildLayer(self,telescope,r0,L0,i_layer,compute_turbulence = True):
         layer.vX            = layer.windSpeed*np.sin(np.deg2rad(layer.direction))      
         
         # Diameter and resolution of the layer including the Field Of View and the number of extra pixels
-        layer.D             = self.telescope.D+2*np.tan(self.telescope.fov/2)*layer.altitude*self.oversampling_factor
+        layer.D             = self.telescope.D+2*np.tan(self.fov_rad/2)*layer.altitude*self.oversampling_factor
         layer.resolution    = int(np.ceil((self.telescope.resolution/self.telescope.D)*layer.D))
         
-        layer.D_fov             = self.telescope.D+2*np.tan(self.telescope.fov/2)*layer.altitude
+        layer.D_fov             = self.telescope.D+2*np.tan(self.fov_rad/2)*layer.altitude
         layer.resolution_fov    = int(np.ceil((self.telescope.resolution/self.telescope.D)*layer.D))
         
         layer.center = layer.resolution//2
@@ -438,7 +441,7 @@ def get_covariance_matrices(self,layer):
         # Compute the covariance matrices
         compute_covariance_matrices = True
 
-        if self.telescope.fov == 0: 
+        if self.fov_rad == 0: 
             try:
                 c=time.time()        
                 self.ZZt_r0 = self.ZZt_r0
@@ -562,14 +565,18 @@ def print_atm(self):
        
     def print_properties(self):
         print('%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ATMOSPHERE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
-        print('{: ^12s}'.format('Layer') + '{: ^12s}'.format('Direction')+ '{: ^12s}'.format('Speed')+ '{: ^12s}'.format('Altitude')+ '{: ^12s}'.format('Cn2') )
-        print('{: ^12s}'.format('') + '{: ^12s}'.format('[deg]')+ '{: ^12s}'.format('[m/s]')+ '{: ^12s}'.format('[m]')+ '{: ^12s}'.format('[m-2/3]') )
+        print('{: ^12s}'.format('Layer') + '{: ^12s}'.format('Direction')+ '{: ^12s}'.format('Speed')+ '{: ^12s}'.format('Altitude')+ '{: ^12s}'.format('Cn2')+ '{: ^12s}'.format('Diameter') )
+        print('{: ^12s}'.format('') + '{: ^12s}'.format('[deg]')+ '{: ^12s}'.format('[m/s]')+ '{: ^12s}'.format('[m]')+ '{: ^12s}'.format('[m-2/3]') + '{: ^12s}'.format('[m]'))
 
-        print('----------------------------------------------------------------------')
+        print('======================================================================')
         
         for i_layer in range(self.nLayer):
-            print('{: ^12s}'.format(str(i_layer+1)) + '{: ^12s}'.format(str(self.windDirection[i_layer]))+ '{: ^12s}'.format(str(self.windSpeed[i_layer]))+ '{: ^12s}'.format(str(self.altitude[i_layer]))+ '{: ^12s}'.format(str(self.fractionalR0[i_layer]) ))
-            print('------------------------------------------------------------------')
+            print('{: ^12s}'.format(str(i_layer+1)) + '{: ^12s}'.format(str(self.windDirection[i_layer]))+ '{: ^12s}'.format(str(self.windSpeed[i_layer]))+ '{: ^12s}'.format(str(self.altitude[i_layer]))+ '{: ^12s}'.format(str(self.fractionalR0[i_layer]) ) + '{: ^12s}'.format(str(getattr(self,'layer_'+str(i_layer+1)).D )))
+            if i_layer<self.nLayer-1:
+                print('----------------------------------------------------------------------')
+
+        print('======================================================================')
+
 
         print('{: ^18s}'.format('r0 @500 nm') + '{: ^18s}'.format(str(self.r0)+' [m]' ))
         print('{: ^18s}'.format('L0') + '{: ^18s}'.format(str(self.L0)+' [m]' ))
@@ -578,8 +585,22 @@ def print_properties(self):
         print('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
         
     def __mul__(self,obj):
-        if obj.tag == 'telescope':
-            self.telescope = obj
+        if obj.tag == 'telescope' or obj.tag == 'source':
+            if obj.tag == 'telescope':
+                if self.fov == obj.fov:   
+                    self.telescope = obj
+                else:          
+                    print('Re-initializing the atmosphere to match the new telescope fov')
+                    self.hasNotBeenInitialized = True
+                    self.initializeAtmosphere(obj)
+
+            else:
+                if obj.coordinates[0] <= self.fov/2:
+                    self.telescope.src = obj
+                    obj = self.telescope
+                else:
+                    raise ValueError('The source object zenith ('+str(obj.coordinates[0])+'") is outside of the telescope fov ('+str(self.fov//2)+'")! You can:\n - Reduce the zenith of the source \n - Re-initialize the atmosphere object using a telescope with a larger fov')
+            
             self.set_pupil_footprint()     
             phase_support = self.initialize_phase_support()
             for i_layer in range(self.nLayer):
@@ -599,7 +620,7 @@ def __mul__(self,obj):
             obj.isPaired     = True
             return obj
         else:
-            raise AttributeError('The atmosphere can be multiplied only with a Telescope object! ')
+            raise AttributeError('The atmosphere can be multiplied only with a Telescope or a Source object!')
     
     
     def display_atm_layers(self,layer_index= None,fig_index = None):
diff --git a/OOPAO/Telescope.py b/OOPAO/Telescope.py
index dce7958..8f8a628 100644
--- a/OOPAO/Telescope.py
+++ b/OOPAO/Telescope.py
@@ -125,7 +125,8 @@ def __init__(self,resolution:float, diameter:float,samplingTime:float=0.001,cent
         self.D                           = diameter                  # Diameter in m
         self.pixelSize                   = self.D/self.resolution    # size of the pixels in m
         self.centralObstruction          = centralObstruction        # central obstruction
-        self.fov                         = fov                       # Field of View 
+        self.fov                         = fov      # Field of View in arcsec converted in radian
+        self.fov_rad                     = fov/206265      # Field of View in arcsec converted in radian
         self.samplingTime                = samplingTime              # AO loop speed
         self.isPetalFree                 = False                     # Flag to remove the petalling effect with ane ELT system. 
         self.index_pixel_petals          = None                      # indexes of the pixels corresponfong to the M1 petals. They need to be set externally
@@ -684,8 +685,13 @@ def print_properties(self):
         print('{: ^18s}'.format('Surface')                      + '{: ^18s}'.format(str(np.round(self.pixelArea*self.pixelSize**2)))    +'{: ^18s}'.format('[m2]'  ))
         print('{: ^18s}'.format('Central Obstruction')          + '{: ^18s}'.format(str(100*self.centralObstruction))                   +'{: ^18s}'.format('[% of diameter]' ))
         print('{: ^18s}'.format('Pixels in the pupil')          + '{: ^18s}'.format(str(self.pixelArea))                                +'{: ^18s}'.format('[pixels]' ))
+        print('{: ^18s}'.format('Field of View')                + '{: ^18s}'.format(str(self.fov))                                      +'{: ^18s}'.format('[arcsec]' ))
         if self.src:
-            print('{: ^18s}'.format('Source '+self.src.type)            + '{: ^18s}'.format(str(np.round(1e9*self.src.wavelength,2)))            +'{: ^18s}'.format('[nm]' ))
+            if self.src.type == 'asterism':
+                for i_src in range(len(self.src.src)):
+                    print('{: ^18s}'.format('Source '+self.src.src[i_src].type)            + '{: ^18s}'.format(str(np.round(1e9*self.src.src[i_src].wavelength,2)))            +'{: ^18s}'.format('[nm]' ))
+            else:
+                print('{: ^18s}'.format('Source '+self.src.type)            + '{: ^18s}'.format(str(np.round(1e9*self.src.wavelength,2)))            +'{: ^18s}'.format('[nm]' ))
         else:
             print('{: ^18s}'.format('Source')            + '{: ^18s}'.format('None')                                             +'{: ^18s}'.format('' ))