FUNCTION SCC_GET_STARS, files, full=full,cor2_point=cor2_point, _extra=_extra ;+ ; $Id: scc_get_stars.pro,v 1.4 2018/08/29 13:17:43 colaninno Exp $ ; ; Project : STEREO SECCHI ; ; Name : scc_get_stars ; ; Purpose : return information for stars in a given secchi image ; ; Explanation: ; ; Use : IDL> result = scc_get_stars(files) ; ; Inputs : string or string array indicating fits file(s) for which star info is desired ; ; Outputs : returns single structure or pointer array to multiple structures ; ; Keywords : ; FULL: Query the full (+35k stars) SECCHI star catalog. Default is the 6k stars ; sub-catalog (where Vmag < 7.0) ; ; Calls from LASCO : NONE ; ; Common : ; ; Restrictions: requires use of IDL SPICE routines. Currently only good for Cor1 or Cor2. ; ; Side effects: Initial reading of _full_ star catalog takes a few seconds ; ; Category : Astrometry (see http://durance.oamp.fr/lasco/fichiers/soft_categories.html) ; ; Prev. Hist. : None. ; ; Written : Karl Battams, NRL/I2, June 2006 ; ; $Log: scc_get_stars.pro,v $ ; Revision 1.4 2018/08/29 13:17:43 colaninno ; changed running cor2_point from default to keyword ; ; Revision 1.3 2008/04/01 17:44:39 reduce ; fixed star search. Now does HI and uses COR2_POINT. Karl B. ; ; Revision 1.2 2008/03/17 20:46:05 avourlid ; Fixed problem for star selection around Sun RA=0 ; ; Revision 1.4 2008/02/15 15:25:05 reduce ; Star pixel positions returned in structure. Karl B. ; ; Revision 1.3 2008/02/06 20:30:00 reduce ; took out 'devpath' kw call to rd_catalog. KarlB ; ; Revision 1.2 2006/09/06 16:04:22 reduce ; Add /full kw for reading full secchi catalog. ; ; Revision 1.1 2006/09/05 20:04:15 reduce ; Added to cvs. Karl B. ; ;- ;================================================================================ ; NOTE: If more than one input file is passed, the value returned will be a ; pointer array to a set of structures. For a procedure call of... ; IDL> stars = scc_get_stars(filelist) ; ... the pointer elements can be accessed as follows: ; IDL> help,/str,(*(stars[0]))[0] ; where the first index is the index of te pointer array and the second ; index is the structure array index. ; Type IDL>?pointers for more information on pointers in IDL. ;================================================================================ IF datatype(files) NE 'STR' THEN BEGIN PRINT,'"Files" parameter must be a string or string array.' return,-1 ENDIF ; Read the appropriate star catalog -- we only want to do this once scc_rd_catalog,catalog,full=full,_extra=_extra nf = N_ELEMENTS(files) IF (nf GT 1) THEN star_ptrarr = ptrarr(nf) FOR i=0,nf-1 DO BEGIN ; Gather parameters needed for finding RA-DEC of Sun hdr = headfits(files[i]) ; just need to read FITS header date = sxpar(hdr,'DATE-OBS') ; ?? ;time = sxpar(hdr,'TIME-OBS') ;dt = strcompress(date,/remove_all)+'T'+strcompress(time,/remove_all) sc = sxpar(hdr,'OBSRVTRY') sc = strcompress(sc,/remove_all) sc = strmid(sc,strlen(sc)-1,1) ;determine which instrument we have: tel = sxpar(hdr,'DETECTOR') tel = strlowcase(strcompress(tel,/remove_all)) IF ((tel NE 'hi1') and (tel NE 'hi2')) THEN BEGIN ; GET SUN RA-DEC for COR2 only sunpos = GET_STEREO_COORD( date, 'SUN', system='gei', /novel ) scpos = GET_STEREO_COORD( date, sc, system='gei', /novel ) state = sunpos-scpos cspice_reclat, state, distance, ra, dec ra = ra * (180.d0 / !dpi) dec = dec * (180.d0 / !dpi) print, 'SUN RA and DEC: ',ra, dec ENDIF ; Next we need to determine the RA-DEC bounds of the fov hdr=scc_fitshdr2struct(hdr) ; Now convert this to a struct for convenience ; Now get RA-DEC _BOX_ bounds (we'll worry about circles later...) ; Note that for HI, we use the center-pixel value as the box center (as opposed to ; Sun-center in COR) CASE strlowcase(tel) of 'cor1' : BEGIN angwidth = 2.0 ; 4.3 deg field for Cor1 END 'cor2' : BEGIN angwidth = 8.6 ; 8.6 deg field for Cor1 IF keyword_set(cor2_point) THEN cor2_point,hdr wcs=fitshead2wcs(hdr,system='A') ; make the wcs header END 'hi1' : BEGIN wcs=fitshead2wcs(hdr,system='A') ; make the wcs header angwidth = 28 ; make this oversized on purpose cenpix=[hdr.naxis1/2,hdr.naxis1/2] cencoord=wcs_get_coord(wcs,cenpix) ra=cencoord[0] dec=cencoord[1] ; Sometimes, for HI1-B, you have to play with the roll angle... ; rot=0.29 ; newang=wcs.roll_angle+rot ; pc=[ [cos(newang/radeg),-sin(newang/radeg)],[sin(newang/radeg),cos(newang/radeg)] ] ; wcs.pc=pc ; wcs.roll_angle=newang END 'hi2' : BEGIN PRINT,'HI2 images not currently supported. Please check back soon...' ; 70-degrees return,-1 END ELSE : BEGIN PRINT,'Camera not found.' return,-1 ENDELSE ENDCASE ; ; Now define the field of view based on sun/image center angrad = 0.5 * angwidth ; angrad is the angular radius of the fov if (ra LT 0) then ra = 360 + ra ; 360 added because spice gives ra in +/- 180 range ; RA increases right-to-left, so we define left_ra as the RA on the left, and right_ra as ; the RA on the right of the image. ; left_ra will USUALLY be greater than right_ra EXCEPT when crossing the 0/360 boundary... left_ra = ra + angrad right_ra = ra - angrad updec = dec + angrad lowdec = dec - angrad ;now get the subset of stars of interest ; simplify things to a box (ignore COR's circular fov for now) IF (left_ra GT right_ra) then BEGIN star_struct = catalog( where( (catalog.alpha LT left_ra) AND (catalog.alpha GT right_ra) AND $ (catalog.delta LT updec) AND (catalog.delta GT lowdec) ) ) ENDIF ELSE BEGIN ; this is if we cross the RA=0 boundary star_struct = catalog( where( (catalog.alpha GT right_ra) OR $ ((catalog.alpha GT 0) AND (catalog.alpha LT left_ra)) AND $ (catalog.delta LT updec) AND (catalog.delta GT lowdec) ) ) ENDELSE ; For COR, apply circular aperture and exclude stars behind occulter: IF ((tel EQ 'cor2') or (tel EQ 'cor1')) THEN BEGIN sunra = ra sundec = dec star_sun_vec=sqrt((star_struct.alpha - sunra)^2 + (star_struct.delta - sundec)^2 ) good = where( (star_sun_vec LE angrad) and (star_sun_vec GT 1.0) ) IF (n_elements(good) EQ 1) THEN BEGIN IF (good EQ -1) THEN BEGIN PRINT,'****** ERROR!!! No stars found in fov! ******** PRINT,'This could be a bug or it could be that there are no *bright* stars PRINT,'in the fov. In that case, you should use the /full keyword to search the PRINT,'full SECCHI catalog. Good luck! :) return,-1 ENDIF ENDIF star_struct=star_struct[good] ENDIF ; COR2_POINT should fix this stuff now so I'll comment it out and remove it eventually (KB) ; goto,skip ;IF ( sc EQ 'A' ) THEN BEGIN ; Calibrate the A-spacecraft images ; radeg=(180/!dpi) ; ;stop ; wcs.crpix=[1017.5,1024.5] ; crd=[0.0,0.0] ; convert_stereo_lonlat,date,crd,'HPC','GEI',spacecraft=sc,/degrees ; wcs.crval=crd ; rot=-0.45 ; newang=wcs.roll_angle+rot ; pc=[ [cos(newang/radeg),-sin(newang/radeg)],[sin(newang/radeg),cos(newang/radeg)] ] ; wcs.pc=pc ; wcs.roll_angle=newang ; wcs.cdelt=[-14.7/3600.0, 14.7/3600.0] ;ENDIF ELSE BEGIN ; Calibrate the B-spacecraft images ; radeg=(180/!dpi) ; ;stop ; wcs.crpix=[1028.7,1025.5] ; crd=[0.0,0.0] ; convert_stereo_lonlat,date,crd,'HPC','GEI',spacecraft=sc,/degrees ; wcs.crval=crd ; rot=0.2 ; newang=wcs.roll_angle+rot ; ;stop ; pc=[ [cos(newang/radeg),-sin(newang/radeg)],[sin(newang/radeg),cos(newang/radeg)] ] ; wcs.pc=pc ; wcs.roll_angle=newang ; wcs.cdelt=[-14.7/3600.0, 14.7/3600.0] ;ENDELSE ;skip: ; END KB testing ;stop ; get pixel positions of the stars aa = [[star_struct.alpha],[star_struct.delta]] aa = transpose(aa) coord = wcs_get_pixel(wcs,aa) nc=n_elements(coord[0,*]) ; ; Add pixels positions to the structure newstr=add_tag(star_struct,0.0,'XPIXEL') newstr=add_tag(newstr,0.0,'YPIXEL') star_struct=newstr FOR k=0,nc-1 DO BEGIN star_struct[k].xpixel=coord[0,k] star_struct[k].ypixel=coord[1,k] ENDFOR ; weed out the stars that are outside the fov still_good = where((star_struct.xpixel GT 2) and (star_struct.xpixel LT (hdr.naxis1 - 2)) and (star_struct.ypixel gt 2) and (star_struct.ypixel LT (hdr.naxis1 - 2) )) star_struct=star_struct[still_good] if nf GT 1 THEN star_ptrarr[i] = ptr_new(star_struct) ENDFOR if nf GT 1 THEN BEGIN PRINT,'*** Returning a pointer array of ',strcompress(nf,/remove_all),' elements. ****' return,star_ptrarr ENDIF ELSE return,star_struct END