COSMIC-PAH / photoelectric-heating

Blame view

coupes_fits.py 8.59 KB
Edit Raw Normal View History
bb489272   Jalabert   initial commit 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
 
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Sun Mar 15 16:37:04 2020
@author: sacha
"""

def cut_fits(filename,coords,plot_all_cut,average,model_output=[],fig_index=1):
    '''Function allowing to make a certain number of cuts between 2 points specified
    - can computes error bars by averaging the surronding pixels
    - plots the map and the corresponding cuts
    
    input:
        filename : str, name of the .fits file in which the cuts are made
        coords : list with the element arranged in 3 parts, where 'n' is the number of cuts:
                    - first list: RA/DEC coordinates of the n starting points, coordinates in deg., each RA/DEC couple is a list,
                    - second list: RA/DEC coordinates of the n ending points, coordinates in deg.,
                    - third item: RA/DEC coordinates of the reference point, ex: a star in the case where the cut starts somewhere else
                    - example with n=3 and the reference being HD200775: coords=[[coord1_1,coord1_2,coord1_3],[coord2_1,coord2_2,coord2_3],coordstar]
        
        plot_all_cut: bool, True if you want to plot the figures,
        average: bool, True if you want the average of a 3x3 square and error bars, False if you just want the pixel values
        model_output: 2xN array, specific to the photoelectric heating effeciency obtained by the model, allowing to plot these values with the cuts.
                      leave it as an empty list if you cut something else than the PE efficiency map.
        fig_index: index of the figure if you do several plots, default to 1
    
    returns list of the data along the cuts, the corresponding errors and the distance in arcsec.
    '''    
    import numpy as np
    import matplotlib.pyplot as plt
    from astropy.io import fits
    from astropy.wcs import WCS
#    from matplotlib.patches import Arrow, Circle    

    # ==============================
    #-- get the data, header and pixel scale in arcsec
    hdu=fits.open(filename)[0]
    data=hdu.data
    hdr=hdu.header
    w_ = WCS(hdr)
    taille_deg0=hdr['CDELT1'] #taille en degres
    taille_sec0=taille_deg0*3600. #bonne taille en arc sec
    
    #-- pixel of reference for the computation of the difference in arcsec
    coordstar=[coords[2][0],coords[2][1]]
    pixs = w_.wcs_world2pix([coordstar],0)
    xs=pixs[0,0].copy()
    ys=pixs[0,1].copy()
    all_data_average=[]
    all_distance_average=[]
    all_error=[]
    cut_list=[]

    #-- copy the data for different purposes
    z = data.copy()
#    zz=z.copy()
    
    if plot_all_cut:
        '''plot of the map to cut'''            
    
        fig = plt.figure(fig_index,figsize=(7.48,5.9))
#        ax = fig.add_subplot(211)
#        pos=ax.imshow(zz[:,::-1],cmap='gray',origin='lower')
#        ax.axis('image')
#        fig.colorbar(pos,ax=ax,extend='both')
        ax2 = fig.add_subplot(111)
        ax2.set_xlabel('distance from the star (")')
        ax2.set_ylabel('$\epsilon_{PAH}$')

    #-- Extract the line...
    # Make a line with "num" points...
    #-- Loop on the number of cuts wanted
    for ind in range(len(coords[0])):
        #-- Extract the coordinate in deg from coords list and convert it into pixel index...
        #... of both extreme points of the cut. From pix0 to pix1
        coord0=[coords[0][ind][0],coords[0][ind][1]]
        coord1=[coords[1][ind][0],coords[1][ind][1]]
        pix0 = w_.wcs_world2pix([coord0],0) # These are in _pixel_ coordinates!!
        pix1 = w_.wcs_world2pix([coord1],0)
        x0=pix0[0,0].copy()
        y0=pix0[0,1].copy()
        x1=pix1[0,0].copy()
        y1=pix1[0,1].copy()


        
        #-- Length of the cut in pixels
        length = int(round(np.sqrt((x0-x1)**2+(y0-y1)**2)))
    #    distance=taille_sec0*np.sqrt((pix_list[:,0]-xs)**2+(pix_list[:,1]-ys)**2)
        
        #-- Creation of vectors x and y of a length of "length" between the extreme positions...
        #... and make sure that they are integers
        x, y = np.linspace(x0, x1, length), np.linspace(y0, y1, length)
        x_,y_=x.round(),y.round()
        x_int,y_int=x_.astype(np.int), y_.astype(np.int)
        #-- coordinates of each pixels
        coord=np.zeros([len(x_int),2])
        coord[:,0]=x_int
        coord[:,1]=y_int
        
        # Extract the values along the line
        distance_average=[]
        data_average=[]
        pixel_value=[]
        error=[]
        list_of_z=[]
        list_of_z2=[]
        plouf=[]
        g=0
        #-- Loop on the pixels inside the cut
        for xi,yi in coord:
            #-- fix the pixel position and extract the value
            xi,yi=xi.astype(np.int), yi.astype(np.int)
            #-- distance from the reference point
            distance=taille_sec0*np.sqrt((xi-xs)**2+(yi-ys)**2)
            distance_average.append(distance)
            
            #-- if you want a value with assymetric errorbars around a mean value
            #... average the value of the pixel with all the surronding pixels
            if average:
              
                if (yi==z.shape[0]-1):
                    list_of_z.append([z[yi,xi],z[yi,xi+1],z[yi-1,xi+1],z[yi-1,xi],z[yi-1,xi-1],z[yi,xi-1]])
                    
                    if np.all(np.isnan(list_of_z)==True):
                        data_average.append(np.mean(np.array([z[yi,xi],z[yi,xi+1],z[yi-1,xi+1],z[yi-1,xi],z[yi-1,xi-1],z[yi,xi-1]])))
                        error.append(np.std(np.array([z[yi,xi],z[yi,xi+1],z[yi-1,xi+1],z[yi-1,xi],z[yi-1,xi-1],z[yi,xi-1]])))
                    else:
                        for plouf in list_of_z[0]:
                            g+=1
                            if np.isnan(plouf)==False:
                                list_of_z2.append(plouf)
    #                            list_of_z2gamma.append(list_of_zgamma[0][g])
                        data_average.append(np.mean(list_of_z2))
    #                    gammaverage.append(np.mean(list_of_z2gamma))
                        error.append(np.std(list_of_z2))                    
                                
                
                else:
                    list_of_z.append([z[yi,xi],z[yi+1,xi-1],z[yi+1,xi],z[yi+1,xi+1],z[yi,xi+1],z[yi-1,xi+1],z[yi-1,xi],z[yi-1,xi-1],z[yi,xi-1]])
    #                list_of_zgamma.append([data_gamma[yi,xi],data_gamma[yi,xi+1],data_gamma[yi-1,xi+1],data_gamma[yi-1,xi],data_gamma[yi-1,xi-1],data_gamma[yi,xi-1]])
    
        
                    if np.all(np.isnan(list_of_z)==True):
                        data_average.append(np.mean(np.array([z[yi,xi],z[yi+1,xi-1],z[yi+1,xi],z[yi+1,xi+1],z[yi,xi+1],z[yi-1,xi+1],z[yi-1,xi],z[yi-1,xi-1],z[yi,xi-1]])))
                        error.append(np.std(np.array([z[yi,xi],z[yi+1,xi-1],z[yi+1,xi],z[yi+1,xi+1],z[yi,xi+1],z[yi-1,xi+1],z[yi-1,xi],z[yi-1,xi-1],z[yi,xi-1]])))
    #                    gammaverage.append(np.mean(np.array([data_gamma[yi,xi],data_gamma[yi,xi+1],data_gamma[yi-1,xi+1],data_gamma[yi-1,xi],data_gamma[yi-1,xi-1],data_gamma[yi,xi-1]])))
    
    
                    else:
                        for plouf in list_of_z[0]:
                            if np.isnan(plouf)==False:
                                list_of_z2.append(plouf)
                        data_average.append(np.mean(list_of_z2))
                        error.append(np.std(list_of_z2))
                        
                
                list_of_z=[]
                list_of_z2=[]
                plouf=[]   
            else:
                pixel_value.append(z[yi,xi])    
                
        #-- save the vector of distances of the current cut
        all_distance_average.append(distance_average)

        if average:
            #-- If you want just the value of the pixel with no error bar
            
            all_data_average.append(np.array(data_average)  )
            all_error.append(error) 


            if plot_all_cut:
                ax2.errorbar(np.array(all_distance_average[ind]),all_data_average[ind],yerr=all_error[ind],label='{}'.format(ind))#,label='cut{}'.format(ind))
                ax2.legend()
        else:
            cut_list.append(pixel_value)
        
          
            if plot_all_cut:
                ax2.semilogy(all_distance_average[ind],cut_list[ind])
    #-- Plot
    if model_output!=[]:
        eff=model_output[1,:].copy()
        dist=model_output[0,:].copy()
        ax2.semilogy(dist,eff,'*',color='red',label='modèle, G=0.5')
        ax2.axvline(42.8,lw=0.8,ls='--',color='blue')
        ax2.axvline(49.7,lw=0.8,ls='--',color='orange')
        ax2.axvline(37.2,lw=0.8,ls='--',color='green')
        ax2.axhline(0.0164,lw=0.8,ls='--',color='green')
        
    plt.legend()


    return all_data_average,all_error,all_distance_average