Changeset 28 for mire


Ignore:
Timestamp:
09/12/08 16:46:03 (12 years ago)
Author:
meynadie
Message:

ajout scripts centrage

Location:
mire
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • mire/mire.sex

    r16 r28  
    1616THRESH_TYPE      RELATIVE       # threshold type: RELATIVE (in sigmas) 
    1717                                # or ABSOLUTE (in ADUs) 
    18 DETECT_THRESH    10            # <sigmas> or <threshold>,<ZP> in mag.arcsec-2 
     18DETECT_THRESH    5             # <sigmas> or <threshold>,<ZP> in mag.arcsec-2 
    1919ANALYSIS_THRESH  10            # <sigmas> or <threshold>,<ZP> in mag.arcsec-2 
    2020  
     
    6161MAG_GAMMA        4.0            # gamma of emulsion (for photographic scans) 
    6262GAIN             0.0            # detector gain in e-/ADU 
    63 PIXEL_SCALE      0.944          # size of pixel in arcsec (0=use FITS WCS info) 
     63PIXEL_SCALE      0              # size of pixel in arcsec (0=use FITS WCS info) 
    6464  
    6565#------------------------- Star/Galaxy Separation ---------------------------- 
  • mire/mire.sh

    r27 r28  
    2121CATALOG=$(basename $1 .fits).cat 
    2222OUTPUT_FITS=$(basename $1 .fits)_sf.fits 
     23OUTPUT_RING=$(basename $1 .fits).ring 
    2324 
    2425if test $skip_sextractor == false 
     
    3334echo "Mesures géométriques..." 
    3435OUTPUT_FILENAME=$(basename $1 .fits).mes 
    35 python $PATH_TO_MIRE/process.py -c $CATALOG -r $2 -o $OUTPUT_FITS -m $OUTPUT_FILENAME -t $3 
     36python $PATH_TO_MIRE/process.py -c $CATALOG -r $2 -o $OUTPUT_FITS -m $OUTPUT_FILENAME -t $3 --ring $OUTPUT_RING 
    3637 
    3738 
  • mire/mire_all.sh

    r27 r28  
    1414 
    1515rm -f collect 
     16rm -f *.mes 
    1617for file in *_0.fits; do 
    1718    if test $skip_sextractor == true 
     
    2425    done; 
    2526 
    26  
    27  
     27moyenne_ring.py > moyenne_ring 
  • mire/process.py

    r27 r28  
    2424                  default="output.mes", 
    2525                  help="Name of the output measurement file") 
     26parser.add_option("--ring",dest="ring", 
     27                  default="output.ring", 
     28                  help="Name of the output measurement file") 
    2629parser.add_option("-t","--timing",dest="timing", 
    2730                  help="Name of the input timing file") 
     
    8487    for source in catalog: 
    8588        pos_source = array([source['X_IMAGE'],source['Y_IMAGE']]) 
    86         if ((modulus(rep[i,:] - pos_source)) < 5): 
     89        max_dist = 5 # pixels 
     90        if ((modulus(rep[i,:] - pos_source)) < max_dist): 
     91            max_dist = modulus(rep[i,:] - pos_source) 
    8792            detections +=1 
    8893            real_rep[i,:] = pos_source 
    89     if (detections !=1): 
     94    if (detections > 1): 
    9095        print "Attention ! ambiguïté pour repÚre ",i,\ 
    9196              ", nb de détections :",detections 
    9297         
     98    if (detections ==0): 
     99        print "Erreur ! Pas de détection pour repÚre ",i 
     100        exit (2) 
    93101 
    94102# Détermination du centre : barycentre des 4 repÚres 
     
    173181            if (nb_voisin > 0): 
    174182                catalog[knot_idx]['MEAN_DIST_TO_NGBRS'] = float(distance)/weights 
    175                 lookup_sources[i,j,3] = float(distance)/weights 
     183                lookup_sources[i,j,3] = 60./ (float(distance)/weights) 
    176184         
    177185cleaned_catalog = [] 
     
    250258 
    251259a = lookup_sources[:,:,3] 
    252              
     260 
    253261hdu = pyfits.PrimaryHDU(a) 
    254262hdulist = pyfits.HDUList([hdu]) 
     
    268276 
    269277meas.write("# Moyenne des mesures /  Module vec_i /  Module vec_j / deltaT (min) / date\n") 
    270 outs = "%7.4f %7.4f %7.4f %7.3f %s\n" % (sum(tmp)/len(tmp), 
    271                                          modulus(vec_i), 
    272                                          modulus(vec_j), 
    273                                          phase, 
    274                                          dt.isoformat()) 
     278# transfo pixels / arcminutes -> mas / pixel 
     279moyenne_mesures = 60 * (1. / (sum(tmp)/len(tmp))) 
     280module_veci = 60 * (1. /  modulus(vec_i)) 
     281module_vecj = 60 * (1. /  modulus(vec_j)) 
     282 
     283outs = "%10.7f %10.7f %10.7f %7.3f %s\n" % (moyenne_mesures, 
     284                                            module_veci, 
     285                                            module_vecj, 
     286                                            phase, 
     287                                            dt.isoformat()) 
    275288meas.write(outs) 
    276289 
    277290meas.close() 
    278291 
     292## Facteur d'échelle par rapport au centre du champ sur une couronne 
     293## représentative du diamÚtre solaire 
     294 
     295# Quel noeud est (à peu prÚs) au centre ? 
     296 
     297 
     298dist_mini = 2048 
     299i_centre = 0 
     300for i in range(len(cleaned_catalog)): 
     301    source = cleaned_catalog[i] 
     302    dist = modulus([source["X_IMAGE"] - 1024, source["Y_IMAGE"] - 1024]) 
     303    if (dist < dist_mini): 
     304        dist_mini = dist 
     305        i_centre = i 
     306 
     307#print "Source au centre : ",cleaned_catalog[i_centre] 
     308 
     309# On parcourt une couronne de noeuds et on calcule le facteur d'échelle pour chacun 
     310 
     311rmin = 14 # en inter-noeud 
     312rmax = 18 # en inter-noeud 
     313 
     314source_centrale = cleaned_catalog[i_centre] 
     315 
     316couronne=[] 
     317 
     318for source in cleaned_catalog: 
     319    tmp = {} 
     320    dist_knots = modulus([source["I_GRID"] - source_centrale["I_GRID"], 
     321                          source["J_GRID"] - source_centrale["J_GRID"]]) 
     322    x_rel = source["X_IMAGE"] - source_centrale["X_IMAGE"] 
     323    y_rel = source["Y_IMAGE"] - source_centrale["Y_IMAGE"] 
     324    dist_pix = modulus([x_rel,y_rel]) 
     325    if ((dist_knots >= rmin) and (dist_knots <= rmax)): 
     326        tmp['ScFact'] = dist_knots * 60. / dist_pix # arcsec / pix 
     327        tmp['r'] = dist_pix 
     328        tmp['theta'] = math.atan(y_rel/x_rel) 
     329        if (x_rel < 0): 
     330            tmp['theta']+=math.pi 
     331        if (tmp['theta'] < 0): 
     332            tmp['theta'] += 2*math.pi 
     333        if (tmp['theta'] >= 2*math.pi): 
     334            tmp['theta'] -= 2*math.pi 
     335        couronne.append(tmp) 
     336 
     337output_couronne = file(options.ring,'w') 
     338 
     339outs = "# phase: %f date: %s\n" % (phase, 
     340                                 dt.isoformat()) 
     341output_couronne.write(outs) 
     342 
     343for source in couronne: 
     344    outs = "%f %f %10.7f\n" % (source['r'], 
     345                               math.degrees(source['theta']), 
     346                               source['ScFact']) 
     347    output_couronne.write(outs) 
     348 
     349 
     350 
     351output_couronne.close() 
     352         
Note: See TracChangeset for help on using the changeset viewer.