Changeset 28 for mire

Timestamp:

09/12/08 16:46:03 (16 years ago)

Author:

meynadie

Message:

ajout scripts centrage

Location:

mire

Files:

• mire.sex (modified) (2 diffs)
• mire.sh (modified) (2 diffs)
• mire_all.sh (modified) (2 diffs)
• process.py (modified) (5 diffs)

mire/mire.sex

@@ -16 +16 @@
 THRESH_TYPE      RELATIVE       # threshold type: RELATIVE (in sigmas)
                                 # or ABSOLUTE (in ADUs)
-DETECT_THRESH    10            # <sigmas> or <threshold>,<ZP> in mag.arcsec-2
+DETECT_THRESH    5             # <sigmas> or <threshold>,<ZP> in mag.arcsec-2
 ANALYSIS_THRESH  10            # <sigmas> or <threshold>,<ZP> in mag.arcsec-2
  
@@ -61 +61 @@
 MAG_GAMMA        4.0            # gamma of emulsion (for photographic scans)
 GAIN             0.0            # detector gain in e-/ADU
-PIXEL_SCALE      0.944          # size of pixel in arcsec (0=use FITS WCS info)
+PIXEL_SCALE      0              # size of pixel in arcsec (0=use FITS WCS info)
  
 #------------------------- Star/Galaxy Separation ----------------------------

mire/mire.sh

@@ -21 +21 @@
 CATALOG=$(basename $1 .fits).cat
 OUTPUT_FITS=$(basename $1 .fits)_sf.fits
+OUTPUT_RING=$(basename $1 .fits).ring
 
 if test $skip_sextractor == false
@@ -33 +34 @@
 echo "Mesures géométriques..."
 OUTPUT_FILENAME=$(basename $1 .fits).mes
-python $PATH_TO_MIRE/process.py -c $CATALOG -r $2 -o $OUTPUT_FITS -m $OUTPUT_FILENAME -t $3
+python $PATH_TO_MIRE/process.py -c $CATALOG -r $2 -o $OUTPUT_FITS -m $OUTPUT_FILENAME -t $3 --ring $OUTPUT_RING
 
 

mire/mire_all.sh

@@ -14 +14 @@
 
 rm -f collect
+rm -f *.mes
 for file in *_0.fits; do
     if test $skip_sextractor == true
@@ -24 +25 @@
     done;
 
-
-
+moyenne_ring.py > moyenne_ring

mire/process.py

@@ -24 +24 @@
                   default="output.mes",
                   help="Name of the output measurement file")
+parser.add_option("--ring",dest="ring",
+                  default="output.ring",
+                  help="Name of the output measurement file")
 parser.add_option("-t","--timing",dest="timing",
                   help="Name of the input timing file")
@@ -84 +87 @@
     for source in catalog:
         pos_source = array([source['X_IMAGE'],source['Y_IMAGE']])
-        if ((modulus(rep[i,:] - pos_source)) < 5):
+        max_dist = 5 # pixels
+        if ((modulus(rep[i,:] - pos_source)) < max_dist):
+            max_dist = modulus(rep[i,:] - pos_source)
             detections +=1
             real_rep[i,:] = pos_source
-    if (detections !=1):
+    if (detections > 1):
         print "Attention ! ambiguïté pour repÚre ",i,\
               ", nb de détections :",detections
         
+    if (detections ==0):
+        print "Erreur ! Pas de détection pour repÚre ",i
+        exit (2)
 
 # Détermination du centre : barycentre des 4 repÚres
@@ -173 +181 @@
             if (nb_voisin > 0):
                 catalog[knot_idx]['MEAN_DIST_TO_NGBRS'] = float(distance)/weights
-                lookup_sources[i,j,3] = float(distance)/weights
+                lookup_sources[i,j,3] = 60./ (float(distance)/weights)
         
 cleaned_catalog = []
@@ -250 +258 @@
 
 a = lookup_sources[:,:,3]
-            
+
 hdu = pyfits.PrimaryHDU(a)
 hdulist = pyfits.HDUList([hdu])
@@ -268 +276 @@
 
 meas.write("# Moyenne des mesures /  Module vec_i /  Module vec_j / deltaT (min) / date\n")
-outs = "%7.4f %7.4f %7.4f %7.3f %s\n" % (sum(tmp)/len(tmp),
-                                         modulus(vec_i),
-                                         modulus(vec_j),
-                                         phase,
-                                         dt.isoformat())
+# transfo pixels / arcminutes -> mas / pixel
+moyenne_mesures = 60 * (1. / (sum(tmp)/len(tmp)))
+module_veci = 60 * (1. /  modulus(vec_i))
+module_vecj = 60 * (1. /  modulus(vec_j))
+
+outs = "%10.7f %10.7f %10.7f %7.3f %s\n" % (moyenne_mesures,
+                                            module_veci,
+                                            module_vecj,
+                                            phase,
+                                            dt.isoformat())
 meas.write(outs)
 
 meas.close()
 
+## Facteur d'échelle par rapport au centre du champ sur une couronne
+## représentative du diamÚtre solaire
+
+# Quel noeud est (à peu prÚs) au centre ?
+
+
+dist_mini = 2048
+i_centre = 0
+for i in range(len(cleaned_catalog)):
+    source = cleaned_catalog[i]
+    dist = modulus([source["X_IMAGE"] - 1024, source["Y_IMAGE"] - 1024])
+    if (dist < dist_mini):
+        dist_mini = dist
+        i_centre = i
+
+#print "Source au centre : ",cleaned_catalog[i_centre]
+
+# On parcourt une couronne de noeuds et on calcule le facteur d'échelle pour chacun
+
+rmin = 14 # en inter-noeud
+rmax = 18 # en inter-noeud
+
+source_centrale = cleaned_catalog[i_centre]
+
+couronne=[]
+
+for source in cleaned_catalog:
+    tmp = {}
+    dist_knots = modulus([source["I_GRID"] - source_centrale["I_GRID"],
+                          source["J_GRID"] - source_centrale["J_GRID"]])
+    x_rel = source["X_IMAGE"] - source_centrale["X_IMAGE"]
+    y_rel = source["Y_IMAGE"] - source_centrale["Y_IMAGE"]
+    dist_pix = modulus([x_rel,y_rel])
+    if ((dist_knots >= rmin) and (dist_knots <= rmax)):
+        tmp['ScFact'] = dist_knots * 60. / dist_pix # arcsec / pix
+        tmp['r'] = dist_pix
+        tmp['theta'] = math.atan(y_rel/x_rel)
+        if (x_rel < 0):
+            tmp['theta']+=math.pi
+        if (tmp['theta'] < 0):
+            tmp['theta'] += 2*math.pi
+        if (tmp['theta'] >= 2*math.pi):
+            tmp['theta'] -= 2*math.pi
+        couronne.append(tmp)
+
+output_couronne = file(options.ring,'w')
+
+outs = "# phase: %f date: %s\n" % (phase,
+                                 dt.isoformat())
+output_couronne.write(outs)
+
+for source in couronne:
+    outs = "%f %f %10.7f\n" % (source['r'],
+                               math.degrees(source['theta']),
+                               source['ScFact'])
+    output_couronne.write(outs)
+
+
+
+output_couronne.close()
+