Skip to content

Commit

Permalink
Added support for occultation chord offsets
Browse files Browse the repository at this point in the history
  • Loading branch information
@matvii
matvii committed Jul 10, 2016
1 parent 350964e commit 0b03a2d
Show file tree
Hide file tree
Showing 17 changed files with 629 additions and 134 deletions.
15 changes: 9 additions & 6 deletions 135_1.ini
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,9 +6,9 @@ SDLevel=2 #Subdivision level, only 0,1,2 currently support
#LMAX=9 #Maximum number of spherical harmonics to use. Number of coefficients is 3(LMAX+1)^2
MinTim=2443846.0 #Zero time, JD
Angles=52,272,8.40060,0 #Rotation angles beta,lambda (using DAMIT convention) and rotation period (in hours) and initial angle.

#FixShape=1
[Optimization]
#UseAOScaling=1
UseAOScaling=1
NumberOfRounds=70 #Number of optimization runs
LCWeight=8 #Weight given to Lightcurve data
AOWeight=2 #Weight given to AO data
Expand All @@ -17,6 +17,7 @@ ConvexWeight=30 #Weight of convex regularization term
AreaWeight=30 #Weight of area regularizion term, with subdiv only
DiAWeight=1 #Weight of dihedral angle regularization, with subdiv only
OctWeight=10 #Octantoid regularization weight, with octantoids only
ChordWeight=1;
Lambda=0.1 #starting lambda in LM optimization

[Data]
Expand All @@ -25,7 +26,7 @@ UseAO=2 #Set to zero if no AO images, otherwise set to n
UseOC=1 #Set to one if occultation is to be used. Zero otherwise

[LC]
LCFile=Hertha/135.rel #File containing LCs, in DAMIT format
LCFile=135.rel #File containing LCs, in DAMIT format

[AO1]
AOFile=Hertha/jittered_1.fits #AO file in fits format. Image size must even, if not, use AOrect and AOsize below
Expand All @@ -37,7 +38,7 @@ PSFRect=51,51
PSFSize=150,150 #PSF must be of same size as the AO image
PixScale=0.009942,0.009942 #Pixel scale. Must be same in both AOFile and PSFFile
SetCamera=Equ #Camera orientation, either ecliptic or equatorial coordinate system
LowFreq=1
#LowFreq=1
#SetCamUp=0,0,1 #Allows fine-tuning of camera orientation, this vector determines camera orientation with the direction vector
[AO2]
AOFile=Hertha/jittered_2.fits #AO file in fits format. Image size must even, if not, use AOrect and AOsize below
Expand All @@ -49,12 +50,14 @@ PSFRect=51,51
PSFSize=150,150 #PSF must be of same size as the AO image
PixScale=0.009942,0.009942 #Pixel scale. Must be same in both AOFile and PSFFile
SetCamera=Equ #Camera orientation, either ecliptic or equatorial coordinate system
LowFreq=1
#LowFreq=1
#SetCamUp=0,0,1 #Allows fine-tuning of camera orientation, this vector determines camera orientation with the direction vector

Weight=2
[OC]
OCFile=Hertha/Hertha.occ #Occultation filename. Uses same format as KOALA
SetCamera=Equ #Camera orientation
#ChordWeightFile=Hertha/ChordWeight.txt
FreeChords=16,17
[Ephm]
EphFile=Hertha/ephm.dat #This file contains the ephemeris information corresponding to observation times in AO images
#Format is LC E0x E0y E0z Ex Ey Ez, where LC is observation time corresponding to AO files (not necessarily in same order
Expand Down
70 changes: 70 additions & 0 deletions 135_oct.ini
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,70 @@
[Shape]
InitEllipsoid=43.56,43.56,35.64 #Initial shape is ellipsoid with semi-axes a,b,c
Nrows=10 #Number of facets is 8nrows^2. For subdiv surfaces, use 4 or 5, for octantoids, 10
#InitShapeFile=shape2.txt #Use initial shape instead. Only with subdivision surfaces
SDLevel=2 #Subdivision level, only 0,1,2 currently supported
LMAX=6 #Maximum number of spherical harmonics to use. Number of coefficients is 3(LMAX+1)^2
MinTim=2443846.0 #Zero time, JD
Angles=52,272,8.40060,0 #Rotation angles beta,lambda (using DAMIT convention) and rotation period (in hours) and initial angle.
#FixShape=1
[Optimization]
UseAOScaling=1
NumberOfRounds=70 #Number of optimization runs
LCWeight=8 #Weight given to Lightcurve data
AOWeight=2 #Weight given to AO data
OCWeight=0.2 #Weight give to OC data
ConvexWeight=30 #Weight of convex regularization term
AreaWeight=30 #Weight of area regularizion term, with subdiv only
DiAWeight=1 #Weight of dihedral angle regularization, with subdiv only
OctWeight=10 #Octantoid regularization weight, with octantoids only
ChordWeight=1;
Lambda=0.1 #starting lambda in LM optimization

[Data]
UseLC=1 #Use lightcurves
UseAO=2 #Set to zero if no AO images, otherwise set to number of AO images
UseOC=1 #Set to one if occultation is to be used. Zero otherwise

[LC]
LCFile=135.rel #File containing LCs, in DAMIT format

[AO1]
AOFile=Hertha/jittered_1.fits #AO file in fits format. Image size must even, if not, use AOrect and AOsize below
AORect=51,51 #Lower left corner of image. Can be used to take subimage of AOFile. If not set, will use whole image
AOSize=150,150 #Subimage size. Start from (1,1), take 100x110 pixel subimage
Date=2454729.060938 #Observation date (JD) without LT correction. If not set, will use MJD-OBS from fitsfile
PSFFile=Hertha/1_Rpsf.fits #PSF file corresponding to AO image. If not set, will assume psf is delta funtion
PSFRect=51,51
PSFSize=150,150 #PSF must be of same size as the AO image
PixScale=0.009942,0.009942 #Pixel scale. Must be same in both AOFile and PSFFile
SetCamera=Equ #Camera orientation, either ecliptic or equatorial coordinate system
#LowFreq=1
#SetCamUp=0,0,1 #Allows fine-tuning of camera orientation, this vector determines camera orientation with the direction vector
[AO2]
AOFile=Hertha/jittered_2.fits #AO file in fits format. Image size must even, if not, use AOrect and AOsize below
AORect=51,51 #Lower left corner of image. Can be used to take subimage of AOFile. If not set, will use whole image
AOSize=150,150 #Subimage size. Start from (1,1), take 100x110 pixel subimage
Date=2456285.770312 #Observation date (JD) without LT correction. If not set, will use MJD-OBS from fitsfile
PSFFile=Hertha/2_Rpsf.fits #PSF file corresponding to AO image. If not set, will assume psf is delta funtion
PSFRect=51,51
PSFSize=150,150 #PSF must be of same size as the AO image
PixScale=0.009942,0.009942 #Pixel scale. Must be same in both AOFile and PSFFile
SetCamera=Equ #Camera orientation, either ecliptic or equatorial coordinate system
#LowFreq=1
#SetCamUp=0,0,1 #Allows fine-tuning of camera orientation, this vector determines camera orientation with the direction vector
Weight=2
[OC]
OCFile=Hertha/Hertha.occ #Occultation filename. Uses same format as KOALA
SetCamera=Equ #Camera orientation
#ChordWeightFile=Hertha/ChordWeight.txt
FreeChords=16,17
[Ephm]
EphFile=Hertha/ephm.dat #This file contains the ephemeris information corresponding to observation times in AO images
#Format is LC E0x E0y E0z Ex Ey Ez, where LC is observation time corresponding to AO files (not necessarily in same order
# as AO1 and A01) E0=[E0x E0y E0z] is the asteroid->Sun vector and E=[Ex Ey Ez] is asteroid->Earth vector

[Output]
ShapeFile=/tmp/outshape.txt #Write the final shape to this file, in the usual format.
AnglesFile=/tmp/param_1 #Writes the final angles to this file.
LCOutputFile=/tmp/lcout.txt #Writes the final lightcurves to this file.

53 changes: 41 additions & 12 deletions Fit_Occ.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,17 +15,20 @@ void AdjFacet(int *tlist,double *vlist,int nfac,int nvert,int *A)
}
}

void Fit_Occ(int *tlist,double *vlist,int nfac,int nvert,double *angles,double *up,double *E,double *V,double *TIME,double *offset,double *chords,int *type,int nchords,double *W,double *dist,double *dx,double *dy,double *dz,double *dangles,double *dtox,double *dtoy)
void Fit_Occ(int *tlist,double *vlist,int nfac,int nvert,double *angles,double *up,double *E,double *V0,double *TIME,double *offset,double *chords,int *type,int nchords,double *W,double *Chordoffset,double *dist,double *dx,double *dy,double *dz,double *dangles,double *dtox,double *dtoy,double* dCOdoff)
{
/*INPUT:
* TIMES nchordsx2 disappearance and appearance times
* offset=[x,y] offset in plane
* chords nchordsx4 coordinates, disappearance and appearance
* Chordoffset nchords vector,containing chord offset in seconds
* OUTPUT:
* dist nchordsx4 matrix, x and y distance ofclosest model points to disappearance and appearance points
* dx,dy,dz 4nchordsxnvert matrix
* dangles 4nchordsx3 matrix
* dtox,dtoy 4*nchordsx1, derivatives wrt offsets */
* dtox,dtoy 4*nchordsx1, derivatives wrt offsets
* dCOdoff 4*nchordsxnchords matrix, derivatives for chord offsets
* */
double time=0;
double M[3][3],MRT[3][3],MRdb[3][3],MRdl[3][3],MRdo[3][3],MRTT[9];
double MRdbT[9],MRdlT[9],MRdoT[9];
Expand All @@ -35,14 +38,18 @@ void Fit_Occ(int *tlist,double *vlist,int nfac,int nvert,double *angles,double *
int ic1,ic2,fa1,fa2;
double clpoint[2],fapoint[2],dclpx[4],dclpy[4],dfapx[4],dfapy[4];
double temp4[4];
double V[3];
double dLx=0,dLy=0; //Chord offset
for(int j=0;j<2*nchords;j++)
time+=TIME[j];
time=time/(2*nchords);
int *Adj=calloc(nvert*nvert,sizeof(int));
AdjFacet(tlist,vlist,nfac,nvert,Adj);
Calculate_Frame_Matrix(E,up,M);
rotate(angles[0],angles[1],angles[2],angles[3],time,R,Rb,Rl,Ro);

//Calculate velocity in camera frame
V[0]=M[0][0]*V0[0]+M[0][1]*V0[1]+M[0][2]*V0[2];
V[1]=M[1][0]*V0[0]+M[1][1]*V0[1]+M[1][2]*V0[2];
transpose(R,Rt); //Transpose, since we rotate the model, not view directions
transpose(Rb,RbT);
transpose(Rl,RlT);
Expand All @@ -65,11 +72,12 @@ zero_array(dz,4*nchords*nvert);
zero_array(dangles,4*nchords*3);
zero_array(dtox,4*nchords);
zero_array(dtox,4*nchords);
zero_array(dCOdoff,4*nchords*nchords);
transpose2(MRT,MRTT);
transpose2(MRdb,MRdbT);
transpose2(MRdl,MRdlT);
transpose2(MRdo,MRdoT);
double *a,*b;
double a[2],b[2];
double *vlist2=calloc(nvert*3,sizeof(double));
double *vlist2b=calloc(nvert*3,sizeof(double));
double *vlist2l=calloc(nvert*3,sizeof(double));
Expand All @@ -82,13 +90,25 @@ double *a,*b;
double el=0;
for(int j=0;j<nchords;j++)
{
if(Chordoffset!=NULL)
{
dLx=V[0]*Chordoffset[j];
dLy=V[1]*Chordoffset[j];
}
else
{
dLx=0;
dLy=0;
}
if(W!=NULL)
w=W[j];
else
w=1;

a=chords+4*j;
b=chords+4*j+2;
a[0]=chords[4*j]+dLx;
a[1]=chords[4*j+1]+dLy;
b[0]=chords[4*j+2]+dLx;
b[1]=chords[4*j+3]+dLy;

find_chord(tlist,vlist2,nfac,nvert,offset,a,b,Adj,cledge,faedge,clpoint,fapoint,&inters,dclpx,dclpy,dfapx,dfapy);
if(type[j]==-1)
{
Expand Down Expand Up @@ -136,10 +156,10 @@ double *a,*b;



dist[4*j]=w*(clpoint[0]-chords[4*j]); //distance in x
dist[4*j+1]=w*(clpoint[1]-chords[4*j+1]); //dist in y
dist[4*j+2]=w*(fapoint[0]-chords[4*j+2]);
dist[4*j+3]=w*(fapoint[1]-chords[4*j+3]);
dist[4*j]=w*(clpoint[0]-chords[4*j]-dLx); //distance in x
dist[4*j+1]=w*(clpoint[1]-chords[4*j+1]-dLy); //dist in y
dist[4*j+2]=w*(fapoint[0]-chords[4*j+2]-dLx);
dist[4*j+3]=w*(fapoint[1]-chords[4*j+3]-dLy);
ic1=cledge[0];
ic2=cledge[1];
fa1=faedge[0];
Expand Down Expand Up @@ -224,7 +244,11 @@ double *a,*b;
el=dfapy[2]*mr13+dfapy[3]*mr23;
el=w*el;
set_el(dz,4*nchords,nvert,el,4*j+3,fa2);

//Chords offset derivatives
set_el(dCOdoff,4*nchords,nchords,-V[0],4*j,j);
set_el(dCOdoff,4*nchords,nchords,-V[1],4*j+1,j);
set_el(dCOdoff,4*nchords,nchords,-V[0],4*j+2,j);
set_el(dCOdoff,4*nchords,nchords,-V[1],4*j+3,j);
//Offset derivatives
dtox[4*j]=w*(dclpx[0]+dclpx[2]);
dtox[4*j+1]=w*(dclpy[0]+dclpy[2]);
Expand Down Expand Up @@ -261,4 +285,9 @@ double *a,*b;
dist[4*j+3]=w*chords[j*4+3];
}
}
free(vlist2);
free(vlist2b);
free(vlist2l);
free(vlist2o);
free(Adj);
}
10 changes: 10 additions & 0 deletions README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,3 +33,13 @@ Viikinkoski, M; Kaasalainen, M.; Durech, J.: *ADAM: a general method for for usi
##Contact
Bug reports, data, feature suggestions and comments are welcome.
Matti Viikinkoski (matti.viikinkoski@gmail.com)

##Updates
#10.7.2016
- Added support for occultation chord offsets
- Added support for optimizing occultation chord offsets

##TODO
- Albedo variation
- Calibrated lightcurves
- Hapke scattering law
19 changes: 13 additions & 6 deletions adam.ini
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,14 +6,17 @@ SDLevel=2 #Subdivision level, only -1,0,1,2 currently
LMAX=5 #Maximum level of spherical harmonics to use. Number of coefficients is 3(LMAX+1)^2. If this is set, Octantoids are used. Typical values 5-9.
MinTim=2449830.78281 #Zero time, JD0
Angles=24,185,5.079177,0 #Rotation angles beta,lambda (using DAMIT convention) and rotation period (in hours) and initial angle
#FixShape=1 #Fix shape parameters, so only offsets, angles etc. are optimized
#FixAngles=1 #Angles are not optimized
[Optimization]
NumberOfRounds=200 #Number of optimization runs
NumberOfRounds=50 #Number of optimization runs
LCWeight=2 #Weight given to Lightcurve data
AOWeight=2 #Weight given to AO data
ConvexWeight=40 #Weight of convex regularization term
AreaWeight=40 #Weight of area regularizion term, with subdiv only
DiAWeight=0.5 #Weight of dihedral angle regularization. Usually 0.3-3.
OctWeight=10 #Octantoid regularization weight, with octantoids only
ChordWeight=1 #Weight of time offset regularization for free chords.
Lambda=1 #starting lambda in LM optimization
#AOOffsetFile=/tmp/AOoffsets #Read initial AO image offsets from a file
#RDOffsetFile=/tmp/RDoffsets #Read initial range-Doppler offset from a file
Expand All @@ -37,7 +40,7 @@ PSFRect=1,1
PSFSize=150,150 #PSF must be of same size as the AO image
PixScale=0.009942,0.009942 #Pixel scale. Must be same in both AOFile and PSFFile
SetCamera=Equ #Camera orientation, either Ecl or Equ. Equ means image is oriented in equatorial coordinate system (default). If not set, uses CD information from fits file (defaults to equ)
LowFreq=1 #Restrict optimization to low frequency part. Useful for images with a little detail.
#LowFreq=1 #Restrict optimization to low frequency part. Useful for images with a little detail.
#SetCamUp=0,0,1 #Allows fine-tuning of camera orientation, direction 'up' of camera
#Weight=0.5 #Additional weighting for images. Weight>1 gives more weight, Weight<1 means less importance.
[AO2]
Expand All @@ -60,17 +63,21 @@ PixScale=0.15,0.5 #Pixel size of image (df,dt), dt in us.
RadarFreq=2300000000 #Radar frequency
LowFreq=1 #Take only low frequency part of image. May lose some detail, but uses only 1/4 of points.
[Ephm]
EphFile=Metis/ephm.dat #This file contains the ephemeris information corresponding to observation times in all images (RD+AO)
EphFile=Metis/ephm.dat #This file contains the ephemeris information corresponding to observation times in all images
#Format is JD E0x E0y E0z Ex Ey Ez, where JD is observation time corresponding to files (not necessarily in same order
# as AO1 and A01) E0=[E0x E0y E0z] is the asteroid->Sun vector and E=[Ex Ey Ez] is asteroid->Earth vector
[OC]
OCFile=Hertha.occ #Occultation filename. Uses same format as KOALA. See README_occultations.txt
OCFile=Hertha.occ #Occultation filename. Uses same format as KOALA
SetCamera=Equ
#OCCOffset=1,2 #Intial OC offset (x,y) in km for each occultation event. Projection origin is moved, not the chords
#ChordWeightFile=/tmp/ChordWeight #File containing weights for individual chords. Space separates values. Do not use commas.
#ChordWeightFile=/tmp/ChordWeight #File containing weights for individual chords.
ChordOffset=0,0.5,0,0 #Offset for individual chords in seconds. First value corresponds to the first chord
FreeChords=2,15,17 #List of chords that are allowed to move. Value 1 corresponds to the first chord in file OCFile.
#DO NOT SET MISSED OBSERVATION CHORD TO FREE! BAD THINGS WILL HAPPEN! ALSO THE LINE DETERMINED BY FREE
#CHORD MUST INTERSECT THE INITIAL SHAPE. OTHERWISE MORE BAD THINGS WILL HAPPEN!
[Output]
ShapeFile=/tmp/OutputShape.txt #Write the final shape to this file, in the usual (DAMIT) format.
AnglesFile=/tmp/Angles.txt #Writes the final angles to this file.
LCOutputFile=/tmp/lcout.txt #Writes the final lightcurves to this file.
ShapeParamFile=/tmp/Pshape.txt #Parameters used for shape representation
ShapeParamFile=/tmp/Pshape.txt #Parameters used for shape
#RDOffsetFile=/tmp/RDoffset.txt #range-Doppler offset output file
Loading

0 comments on commit 0b03a2d

Please sign in to comment.