Commit 912a3851 authored by Xin Zhang's avatar Xin Zhang
Browse files

Initial commit

parents

.DS_Store

0 → 100644
+12 KiB

File added.

No diff preview for this file type.

View file

getInputPointing.py

0 → 100755
+200 −0
Original line number Diff line number Diff line 
from tkinter.tix import INTEGER

import astropy.coordinates as coord

from astropy import units as u

from pylab import *

import numpy as np

import galsim

# from numba import jit



class Chip(object):

    def __init__(self, chipID):

        self.chipID = chipID



        self.nchip_x = 6

        self.nchip_y = 5

        self.npix_tot_x = 59516

        self.npix_tot_y = 49752

        self.npix_gap_x = (534, 1309)

        self.npix_gap_y = 898



        self.cen_pix_x = 0

        self.cen_pix_y = 0



        self.npix_x = 9216

        self.npix_y = 9232

        self.pix_scale  = 0.074



    def getTanWCS(self, ra, dec, img_rot, pix_scale = None, xcen=None, ycen=None, logger=None):

        """ Get the WCS of the image mosaic using Gnomonic/TAN projection



        Parameter:

            ra, dec:    float

                        (RA, Dec) of pointing of optical axis

            img_rot:    galsim Angle object

                        Rotation of image

            pix_scale:  float

                        Pixel size in unit of as/pix

        Returns:

            WCS of the focal plane

        """

        if logger is not None:

            logger.info("    Construct the wcs of the entire image mosaic using Gnomonic/TAN projection")

        if (xcen == None) or (ycen == None):

            xcen = self.cen_pix_x

            ycen = self.cen_pix_y

        if pix_scale == None:

            pix_scale = self.pix_scale

        dudx =  -np.cos(img_rot.rad) * pix_scale

        dudy =  -np.sin(img_rot.rad) * pix_scale

        dvdx =  -np.sin(img_rot.rad) * pix_scale

        dvdy =  +np.cos(img_rot.rad) * pix_scale

        

        # dudx =  +np.sin(img_rot.rad) * pix_scale

        # dudy =  +np.cos(img_rot.rad) * pix_scale

        # dvdx =  -np.cos(img_rot.rad) * pix_scale

        # dvdy =  +np.sin(img_rot.rad) * pix_scale

        moscen = galsim.PositionD(x=xcen, y=ycen)

        sky_center = galsim.CelestialCoord(ra=ra*galsim.degrees, dec=dec*galsim.degrees)

        affine = galsim.AffineTransform(dudx, dudy, dvdx, dvdy, origin=moscen)

        WCS = galsim.TanWCS(affine, sky_center, units=galsim.arcsec)



        return WCS

    

    def getChipRowCol(self, chipID):

        rowID = ((chipID - 1) % 5) + 1

        colID = 6 - ((chipID - 1) // 5)

        return rowID, colID

    

    def getChipCenter(self):

        """Calculate the edges in pixel for a given CCD chip on the focal plane

        NOTE: There are 5*4 CCD chips in the focus plane for photometric observation.

        Parameters:

            chipID:         int

                            the index of the chip

        Returns:

            A galsim BoundsD object

        """







        chipID = self.chipID



        rowID, colID = self.getChipRowCol(chipID)

        gx1, gx2 = self.npix_gap_x

        gy = self.npix_gap_y



        # xlim of a given CCD chip

        xrem = 2*(colID - 1) - (self.nchip_x - 1)

        xcen = (self.npix_x//2 + gx1//2) * xrem

        if chipID >= 26 or chipID == 21:

            xcen = (self.npix_x//2 + gx1//2) * xrem - (gx2-gx1)

        if chipID <= 5 or chipID == 10:

            xcen = (self.npix_x//2 + gx1//2) * xrem + (gx2-gx1)

        # nx0 = xcen - self.npix_x//2 + 1

        # nx1 = xcen + self.npix_x//2



        # ylim of a given CCD chip

        yrem = (rowID - 1) - self.nchip_y // 2

        ycen = (self.npix_y + gy) * yrem

        # ny0 = ycen - self.npix_y//2 + 1

        # ny1 = ycen + self.npix_y//2



        return galsim.PositionD(xcen, ycen)









# @jit()

def getSelectPointingList(center = [60,-40], radius = 2):

    points = np.loadtxt('skyMapOrSurveyList/sky.dat')



    

    center = center#ra dec

    radius = radius # degree

    c_eclip = coord.SkyCoord(points[:,2]*u.degree, points[:,1]*u.degree,frame='barycentrictrueecliptic')

    c_equtor = c_eclip.transform_to('icrs')



    ids1 = (c_equtor.ra*u.degree).value > center[0]-radius

    ids2 = (c_equtor[ids1].ra*u.degree).value < center[0]+radius

    ids3 = (c_equtor[ids1][ids2].dec*u.degree).value > center[1]-radius

    ids4 = (c_equtor[ids1][ids2][ids3].dec*u.degree).value < center[1]+radius



    num = points[ids1][ids2][ids3][ids4].shape[0]



    p_result = np.zeros([num, 5])

    i = 0



    for p,p_ in zip(points[ids1][ids2][ids3][ids4],c_equtor[ids1][ids2][ids3][ids4]):

        ra = (p_.ra*u.degree).value

        dec = (p_.dec*u.degree).value

        # print(ra, dec)

        lon = p[2]

        lat = p[1]



        p_result[i,0] = ra

        p_result[i,1] = dec

        p_result[i,2] = lon

        p_result[i,3] = lat

        p_result[i,4] = -113.4333

        i = i + 1

    

    return p_result





def findPointingbyChipID(chipID = 8, ra = 60., dec = -40.):

    chip_center = [ra, dec]

    p_list = getSelectPointingList(center = chip_center)

    pchip = Chip(chipID)



    p_num = p_list.shape[0]



    min_d = 1000000000



    r_ra = ra

    r_dec = dec

    for i in np.arange(0,p_num,1):

        ra_n = p_list[i,0]

        dec_n = p_list[i,1]

        rot = p_list[i,4]*galsim.degrees

        chip_wcs = pchip.getTanWCS(ra_n, dec_n, rot)



        c_center = pchip.getChipCenter()



        c_world = chip_wcs.toWorld(c_center)



        ra_s = c_world.ra.deg

        dec_s = c_world.dec.deg

        # print(ra, dec, ra_s, dec_s)

        d = (ra_s - ra)*(ra_s - ra) + (dec_s - dec)*(dec_s - dec)

        if d < min_d:

            min_d = d

            r_ra = ra_n

            r_dec = dec_n



    return galsim.CelestialCoord(ra=r_ra*galsim.degrees,dec=r_dec*galsim.degrees)











if __name__ == "__main__":



    pointing = findPointingbyChipID()



    # chipID = 8

    # pchip = Chip(chipID)

    # ra = pointing.ra.deg

    # dec = pointing.dec.deg

    # rot = -113.433

    # chip_wcs = pchip.getTanWCS(ra, dec, rot*galsim.degrees)



    # c_center = pchip.getChipCenter()



    # c_world = chip_wcs.toWorld(c_center)



    # ra_s = c_world.ra.deg

    # dec_s = c_world.dec.deg

getPointingList.py

0 → 100755
+184 −0
Original line number Diff line number Diff line 
import astropy.coordinates as coord

from astropy import units as u

from pylab import *

import numpy as np

# from numba import jit



def loadSatOrbitDat(datDir='',fileNum=50):

    oData = loadtxt(datDir+'1.txt')



    for i in arange(2, fileNum + 1, 1):

        tdat = loadtxt(datDir+str(i)+'.txt')

        oData = np.vstack((oData,tdat))

    

    return oData







def locateSat(time=2459766., OrbitData = None, startId = 0, orbDataLen = 10000):



    satPos = np.zeros(3)

    nextSid = startId

    for i in np.arange(startId, orbDataLen-1, 1):

        t1 = OrbitData[i,0]

        t2 =  OrbitData[i+1,0]

        if time>= t1 and time<t2:

            if((t2-t1)>130.0/86400.0): # should be 120s, for error, set 120+10s

                     break



            x1 = OrbitData[i,1];

            y1 = OrbitData[i,2];

            z1 = OrbitData[i,3];

            x2 = OrbitData[i+1,1];

            y2 = OrbitData[i+1,2];

            z2 = OrbitData[i+1,3];

			

            l1 = np.sqrt(x1*x1+y1*y1+z1*z1);

            l2 = np.sqrt(x2*x2+y2*y2+z2*z2);

            theta = np.arccos((x1*x2+y1*y2+z1*z2)/(l1*l2));

            theta1 = (time-t1)/(t2-t1)*theta;

            theta2 = theta-theta1;

            l = (t2-time)/(t2-t1)*l1+(time-t1)/(t2-t1)*l2;

			# // double ef = sin(theta2)/sin(theta1);

            

            x0 = np.sin(theta2)*x1/l1+np.sin(theta1)*x2/l2;

            y0 = np.sin(theta2)*y1/l1+np.sin(theta1)*y2/l2;

            z0 = np.sin(theta2)*z1/l1+np.sin(theta1)*z2/l2;

            l_ = np.sqrt(x0*x0+y0*y0+z0*z0);

            

            satPos[0] = x0*l/l_;

            satPos[1] = y0*l/l_;

            satPos[2] = z0*l/l_;

            nextSid = i

            break;	

    return satPos, nextSid









# @jit()

def producePointingList(out = '' , center = [60,-40], radius = 5):

    points = np.loadtxt('skyMapOrSurveyList/sky.dat')



    num = points.shape[0]

    center = center#ra dec

    radius = radius # degree

    c_eclip = coord.SkyCoord(points[:,2]*u.degree, points[:,1]*u.degree,frame='barycentrictrueecliptic')

    c_equtor = c_eclip.transform_to('icrs')



    ids1 = (c_equtor.ra*u.degree).value > center[0]-radius

    ids2 = (c_equtor[ids1].ra*u.degree).value < center[0]+radius

    ids3 = (c_equtor[ids1][ids2].dec*u.degree).value > center[1]-radius

    ids4 = (c_equtor[ids1][ids2][ids3].dec*u.degree).value < center[1]+radius



    elip_txt = '# ra    dec    lon    lat    angle \n'



    for p,p_ in zip(points[ids1][ids2][ids3][ids4],c_equtor[ids1][ids2][ids3][ids4]):

        ra = (p_.ra*u.degree).value

        dec = (p_.dec*u.degree).value

        # print(ra, dec)

        lon = p[2]

        lat = p[1]

        elip_txt = elip_txt + str(round(ra,6)) + '    ' + str(round(dec,6)) + '    ' + str(round(lon,6)) + '    ' + str(round(lat,6)) + '    -113.4333'

        elip_txt = elip_txt + '\n'



    # for i in np.arange(0,num,1):

    #     # lon,lat

    

    #     ra = (c_equtor[i].ra*u.degree).value

    #     dec = (c_equtor[i].dec*u.degree).value



    #     if center[0]-radius <ra < center[0]+radius and center[1]-radius<dec<center[1]+radius:

    #         elip_txt = elip_txt + str(round(ra,6)) + '    ' + str(round(dec,6)) + '    ' + str(round(points[i,2],6)) + '    ' + str(round(points[i,1],6)) + '    -113.4333\n' 



    pointfn = open(out,'w')

    pointfn.write(elip_txt)

    pointfn.flush()

    pointfn.close()





def producePointingList2(out = '' , center = [60,-40], radius = 5, survey_file = 'skyMapOrSurveyList/E17.5_b17.5_beta_11.6_opt_transtime_1_CMG_1_dp_2_0.25_da_10_Texp_1.5_DEC60_500_0.1_800_1000_+5deg.dat'):

    points = np.loadtxt(survey_file)



    num = points.shape[0]

    center = center#ra dec

    radius = radius # degree

    c_eclip = coord.SkyCoord(points[:,2]*u.degree, points[:,1]*u.degree,frame='barycentrictrueecliptic')

    c_equtor = c_eclip.transform_to('icrs')

    ids1 = (c_equtor.ra*u.degree).value > center[0]-radius

    ids2 = (c_equtor[ids1].ra*u.degree).value < center[0]+radius

    ids3 = (c_equtor[ids1][ids2].dec*u.degree).value > center[1]-radius

    ids4 = (c_equtor[ids1][ids2][ids3].dec*u.degree).value < center[1]+radius



    obDataDir = 'orbit20160925/'

    orbitDat = loadSatOrbitDat(obDataDir,50)

    sOrbitId = 0



    elip_txt = '# ra    dec    lon(ecliptic)    lat(ecliptic)    pos_angle    time(julian)    sat_x    sat_y    sat_z    sun_x    sun_y    sun_z    moon_x    moon_y    moon_z    sat_vx    sat_vy    sat_vz    exp_time    isDeep\n'



    for p,p_ in zip(points[ids1][ids2][ids3][ids4],c_equtor[ids1][ids2][ids3][ids4]):

        ra = (p_.ra*u.degree).value

        dec = (p_.dec*u.degree).value

        # print(ra, dec)

        lon = p[2]

        lat = p[1]

        elip_txt = elip_txt + str(round(ra,6)) + '    ' + str(round(dec,6)) + '    ' + str(round(lon,6)) + '    ' + str(round(lat,6)) + '    -113.4333'



        oTime = p[0]

        satPos, sOrbitId =locateSat(time=oTime, OrbitData = orbitDat, startId = sOrbitId, orbDataLen = orbitDat.shape[0])

        tempOrbitId = sOrbitId

        deltT = 0.1 # unit s

        nTime = p[0]+deltT/86400.

        satPosN, _ =locateSat(time=nTime, OrbitData = orbitDat, startId = sOrbitId, orbDataLen = orbitDat.shape[0])

        sat_v = (satPosN - satPos)/deltT

        mm1 = np.sqrt(np.sum(sat_v*sat_v))

        mm2 = np.sqrt(np.sum(satPos*satPos))

        # print(np.sum(sat_v*satPos)/(mm1*mm2))

        # print(sat_v,p[6:9])

        if(np.abs(np.sum(sat_v*satPos)/(mm1*mm2))>0.01):

            print(oTime)

            break





        

        # checkMod = np.sqrt((satPos[0]-p[3])*(satPos[0]-p[3]) + (satPos[1]-p[4])*(satPos[1]-p[4]) + (satPos[2]-p[5])* (satPos[2]-p[5]))



        # if (checkMod > 0.01):

        #     print(oTime, checkMod)

        #     break



        # print(satPos[0]-p[3], satPos[1]-p[4], satPos[2]-p[5])



        # if oTime > 2459768:

        #     break

        

        dat_col = [0,3,4,5,6,7,8,9,10,11,16,14]

        for col in dat_col[:-2]:

            elip_txt = elip_txt + '    ' + str(p[col])

        elip_txt = elip_txt + '    ' + str(sat_v[0]) + '    ' + str(sat_v[1]) + '    ' + str(sat_v[2])



        for col in dat_col[-2:]:

            elip_txt = elip_txt + '    ' + str(p[col])





        elip_txt = elip_txt + '\n'

        



    pointfn = open(out,'w')

    pointfn.write(elip_txt)

    pointfn.flush()

    pointfn.close()





if __name__ == "__main__":

    isRealSurvey = False

    survey_file = 'skyMapOrSurveyList/E17.5_b17.5_beta_11.6_opt_transtime_1_CMG_1_dp_2_0.25_da_10_Texp_1.5_DEC60_500_0.1_800_1000_+5deg.dat'

    outFileName = 'pointing_test_false.dat'

    radi = 1



    center_pos = [60, -40]

    if isRealSurvey:

        producePointingList2(out = outFileName,center = center_pos, radius = radi, survey_file = survey_file)

    else:

        producePointingList(out = outFileName,center = center_pos, radius = radi)

orbit20160925/1.txt

0 → 100755
+30490 −0

File added.

Preview size limit exceeded, changes collapsed.

orbit20160925/10.txt

0 → 100755
+57135 −0

File added.

Preview size limit exceeded, changes collapsed.