LAL: InspiralSpinBankwNDTemplateBank.c Source File

00001 /*
00002 *  Copyright (C) 2007 Chad Hanna, Benjamin Owen
00003 *
00004 *  This program is free software; you can redistribute it and/or modify
00005 *  it under the terms of the GNU General Public License as published by
00006 *  the Free Software Foundation; either version 2 of the License, or
00007 *  (at your option) any later version.
00008 *
00009 *  This program is distributed in the hope that it will be useful,
00010 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
00011 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00012 *  GNU General Public License for more details.
00013 *
00014 *  You should have received a copy of the GNU General Public License
00015 *  along with with program; see the file COPYING. If not, write to the
00016 *  Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston,
00017 *  MA  02111-1307  USA
00018 */
00019 
00020 /**** <lalVerbatim file="InspiralSpinBankwNDTemplateBankCV">
00021  * Authors: Hanna, C. R. and Owen, B. J.
00022  * $Id: InspiralSpinBankwNDTemplateBank.c,v 1.3 2007/06/08 14:41:42 bema Exp $
00023  **** </lalVerbatim> */
00024 
00025 /**** <lalLaTeX>
00026  *
00027  * \subsection{Module \texttt{InspiralSpinBankwNDTemplateBank.c}}
00028  *
00029  * This module creates a bank of templates to search for precessing
00030  * binaries.
00031  *
00032  * \subsubsection*{Prototypes}
00033  * \input{InspiralSpinBankwNDTemplateBankCP}
00034  * %% \idx{LALInspiralSpinBankwNDTemplateBank()}
00035  *
00036  * \subsubsection*{Description}
00037  *
00038  * This function creates a bank of templates to search for precessing
00039  * binaries. It uses the general tiling algorithm LALNDTemplateBank()
00040  *
00041  * \subsubsection*{Algorithm}
00042  *
00043  * The target region of parameter space is a distorted box in the
00044  * coordinates $(x=\psi_0, y=\psi_3, z=\beta)$. The metric at high values of
00045  * $\beta$ is flat. It is convenient to rotate to coordinates $(x',y',z')$
00046  * which lie along eigenvectors of the metric.
00047  *
00048  * The algorithm first draws a rectilinear box in the primed coordinates
00049  * which includes the distorted box, then steps through along the directions
00050  * of the primed coordinates.  At each point it tests if the point lies
00051  * within the distorted box. If the point is inside the distorted box, the
00052  * algorithm adds a template to the linked list. If not, it continues.
00053  *
00054  *
00055  * At the end it copies the linked list into the inspiral package's array
00056  * format.
00057  *
00058  * \subsubsection*{Uses}
00059  *
00060  * \begin{verbatim}
00061  * LALCalloc()
00062  * LALFree()
00063  * \end{verbatim}
00064  *
00065  * \subsubsection*{Notes}
00066  *
00067  *
00068  * The metric relies on approximations that make it valid only for a binary
00069  * system with a total mass $<15M\odot$ where the larger body's minimum mass
00070  * is at least twice the smaller body's maximum mass.  Using 
00071  * that violate these conditions will result in an error message.   
00072  *
00073  * The issue of linked lists vs.\ arrays needs to be seriously addressed. As
00074  * our experience with this code shows, multidimensional tiling of
00075  * complicated parameter regions demands the flexibility of linked lists.
00076  *
00077  * \vfill{\footnotesize\input{InspiralSpinBankwNDTemplateBankCV}}
00078  *
00079  **** </lalLaTeX> */
00080 
00081 
00082 #include <math.h>
00083 #include <lal/AVFactories.h>
00084 #include <lal/FlatMesh.h>
00085 #include <lal/LALConfig.h>
00086 #include <lal/LALConstants.h>
00087 #include <lal/LALDatatypes.h>
00088 #include <lal/LALInspiralBank.h>
00089 #include <lal/LALMalloc.h>
00090 #include <lal/LALStatusMacros.h>
00091 #include <lal/LALStdlib.h>
00092 #include <lal/LIGOMetadataTables.h>
00093 #include <lal/MatrixUtils.h>
00094 #include <lal/SeqFactories.h>
00095 #include <lal/TemplateBankGeneration.h>
00096 
00097 
00098 NRCSID(INSPIRALSPINBANKWNDTEMPLATEBANKC, "$Id: InspiralSpinBankwNDTemplateBank.c,v 1.3 2007/06/08 14:41:42 bema Exp $");
00099 
00100 
00101 /* LALINSPIRALSPINBANKMETRIC() --------------------------------------------- */
00102 void
00103 LALInspiralSpinBankMetric(
00104    LALStatus            *status,
00105    NDTemplateBankInput  *input,
00106    REAL4Array           *metric
00107    )
00108 {
00109   InspiralMomentsEtc moments;           /* Added for LALGetInspiralMoments() */
00110   InspiralTemplate inspiralTemplate;    /* Added for LALGetInspiralMoments() */
00111 
00112   INT2 loop = 0;
00113   REAL4 f0 = input->f0;
00114   REAL8 J1  = 0.0;
00115   REAL8 J4  = 0.0;
00116   REAL8 J6  = 0.0;
00117   REAL8 J9  = 0.0;
00118   REAL8 J11 = 0.0;
00119   REAL8 J12 = 0.0;
00120   REAL8 J14 = 0.0;
00121   REAL8 J17 = 0.0;
00122 
00123   INITSTATUS( status, "LALInspiralSpinBankMetric", INSPIRALSPINBANKWNDTEMPLATEBANKC );
00124   ATTATCHSTATUSPTR( status );
00125   
00126   
00127   inspiralTemplate.fLower  = 30;        /* These are arbitrarily chosen for now */
00128   inspiralTemplate.fCutoff = 2000;      /* They are necessary for LALInspiralGetMoments() */
00129   
00130   LALGetInspiralMoments( status->statusPtr, &moments, input->PSD, &inspiralTemplate );
00131 
00132   /* Rescale the moments to F0 = Noise Curve Minimum */
00133   for(loop = 1; loop <=17; loop++){
00134     moments.j[loop] *= pow((inspiralTemplate.fLower/(f0)), ((7.0-(REAL4) loop)/3.0));
00135     }
00136 
00137 /* This just copies the noise moment data from *moments */
00138   J1  = moments.j[1];
00139   J4  = moments.j[4];
00140   J6  = moments.j[6];
00141   J9  = moments.j[9];
00142   J11 = moments.j[11];
00143   J12 = moments.j[12];
00144   J14 = moments.j[14];
00145   J17 = moments.j[17];
00146                                                                                                                                                 
00147   /* Set metric components as functions of moments. */
00148   metric->data[0] = (REAL4) (1.5)*(J17-J12*J12-(J9-J4*J12)*(J9-J4*J12)/(J1-J4*J4));
00149   metric->data[1] = (REAL4) (1.5)*(J14-J9*J12-(J6-J4*J9)*(J9-J4*J12)/(J1-J4*J4));
00150   metric->data[2] = (REAL4) 0.0;
00151   metric->data[3] = (REAL4) (1.5)*(J14-J9*J12-(J6-J4*J9)*(J9-J4*J12)/(J1-J4*J4));
00152   metric->data[4] = (REAL4) (1.5)*(J11-J9*J9-(J6-J4*J9)*(J6-J4*J9)/(J1-J4*J4));
00153   metric->data[5] = (REAL4) 0.0;
00154   metric->data[6] = (REAL4) 0.0;
00155   metric->data[7] = (REAL4) 0.0;
00156   metric->data[8] = (REAL4) J11-J9*J9-(J6-J4*J9)*(J6-J4*J9)/(J1-J4*J4);
00157   
00158   DETATCHSTATUSPTR( status );
00159   RETURN( status );
00160 } /* LALInspiralSpinBankMetric */
00161 
00162 
00163 /* LALInspiralSpinBankwNDTemplateBank() --------------------------------------------------- */
00164 /* <lalVerbatim file="LALInspiralSpinBankwNDTemplateBankCP"> */
00165 void
00166 LALInspiralSpinBankwNDTemplateBank(
00167     LALStatus            *status,
00168     InspiralTemplateList **tiles,
00169     INT4                 *ntiles,
00170     InspiralCoarseBankIn  coarseIn
00171     )
00172 /* </lalVerbatim> */
00173 {
00174   
00175   /* Initialize variables */
00176   REAL4Array *metric =            NULL; /* parameter-space metric */
00177   UINT4Vector *metricDimensions = NULL; /* contains the dimension of metric */
00178   NDTemplateBankInput NDinput;
00179   NDTemplateBankOutput *NDoutput = NULL;
00180   NDTemplateBankOutput *NDFirst  = NULL;
00181   NDTemplateBankFunctionPtrs NDFunctionPtrs;
00182   INT4 cnt = 0;
00183   REAL4 f0, m1Max, m1Min, m2Max, m2Min;
00184  
00185   /* Set up status pointer. */
00186   INITSTATUS( status, "LALInspiralSpinBankwNDTemplateBank", INSPIRALSPINBANKWNDTEMPLATEBANKC );
00187   ATTATCHSTATUSPTR( status );
00188   
00189   /* Check to make sure that all the parameters are okay */
00190   if (coarseIn.mmCoarse <= 0){
00191     ABORT(status, LALINSPIRALBANKH_ECHOICE, LALINSPIRALBANKH_MSGECHOICE);
00192     }
00193 
00194   if ((coarseIn.mMin <= 0) || (coarseIn.MMax <= 0) || 
00195       (coarseIn.mMin >= coarseIn.MMax) || (3.0*coarseIn.MMax >= 15.0)){
00196     ABORT(status, LALINSPIRALBANKH_ECHOICE, LALINSPIRALBANKH_MSGECHOICE);
00197     }
00198 
00199   /*These parameters have not been added to InspiralCoarseBankIn yet, but when they are the will need to be checked */
00200   /*
00201     if (coarseIn.betaMax < 0) 
00202       ABORT(status, LALINSPIRALBANKH_ECHOICE, LALINSPIRALBANKH_MSGECHOICE);
00203   */
00204 
00205   /* Get 3x3 parameter-space metric. */
00206   /* BEN: mess creating all these structures & adding TRYs etc */
00207   /* BEN: do it by hand, since it's so simple? */
00208 
00209   /* allocate memory for the metric */
00210   LALU4CreateVector( status->statusPtr, &metricDimensions, (UINT4) 2 ); 
00211   metricDimensions->data[0] = 3;
00212   metricDimensions->data[1] = 3;
00213   LALSCreateArray( status->statusPtr, &metric, metricDimensions );
00214 
00215 
00216   /* Hardcode mass range etc for the moment. */
00217   m2Min = NDinput.minParameters[1] = coarseIn.mMin*LAL_MTSUN_SI;
00218   m2Max = NDinput.maxParameters[1] = coarseIn.MMax*LAL_MTSUN_SI;
00219   m1Min = NDinput.minParameters[0] = 2.0*m2Max;
00220   m1Max = NDinput.maxParameters[0] = 15.0*LAL_MTSUN_SI - m2Max;
00221   f0 = NDinput.f0 = 153.0; /*FIX THIS FIX THIS FIX THIS FIX THIS FIX THIS */
00222   
00223   /* Set NDinput parameters */
00224   NDinput.mm = coarseIn.mmCoarse;
00225   NDinput.type = PrecessingType;
00226   NDinput.dimension = 3;
00227   NDinput.PSD = &coarseIn.shf;
00228 
00229   /* Set box on unprimed coordinates including region. */
00230   /*psi0*/ NDinput.minCoordinates[0] = 0.9*(3.0/128) / (pow(LAL_PI*f0*(m1Max+m2Max),1.666667)*(m1Max*m2Max/pow(m1Max+m2Max,2)));
00231   /*psi3*/ NDinput.minCoordinates[1] = 1.1*(-.375*LAL_PI) / (pow(LAL_PI*f0*(m1Max+m2Min),0.6666667)*(m1Max*m2Min/pow(m1Max+m2Min,2)));
00232   /*beta*/ NDinput.minCoordinates[2] = 0;
00233   /*psi0*/ NDinput.maxCoordinates[0] = 1.1*(3.0/128) / (pow(LAL_PI*f0*(m1Min+m2Min),1.666667)*(m1Min*m2Min/pow(m1Min+m2Min,2)));
00234   /*psi3*/ NDinput.maxCoordinates[1] = .9*(-.375*LAL_PI) / (pow(LAL_PI*f0*(m1Min+m2Max),0.6666667)*(m1Min*m2Max/pow(m1Min+m2Max,2)));
00235   /*beta*/ NDinput.maxCoordinates[2] = 3.8* LAL_PI/29.961432 * (1+0.75*m2Max/m1Min) * (m1Max/m2Min) * pow(LAL_MTSUN_SI*100.0/(m1Min+m2Min), 0.6666667);
00236 
00237   /* Set the function pointers for NDTemplate */
00238   NDFunctionPtrs.metric = LALInspiralSpinBankMetric;
00239   NDFunctionPtrs.test = LALInspiralSpinBankBoundary;
00240   
00241   printf("\ncalling LALNDTemplateBank()...\n");
00242   LALNDTemplateBank(status->statusPtr, &NDinput, &NDFunctionPtrs, &NDFirst);
00243   
00244   printf("\nconverting output...\n");
00245   NDoutput = NDFirst;
00246   while(NDoutput->next){
00247     (*ntiles)++;
00248     NDoutput = NDoutput->next;
00249     }
00250 
00251   printf("\n...counter is %i\n", *ntiles);
00252 
00253   *tiles = (InspiralTemplateList *) LALCalloc( *ntiles, sizeof(InspiralTemplateList));
00254   NDoutput = NDFirst;
00255   cnt = 0;
00256   
00257   for (cnt = 0; cnt < *ntiles; cnt++)
00258   {
00259     (*tiles)[cnt].params.mass1 = NDoutput->parameterVals[0];
00260     (*tiles)[cnt].params.mass2 = NDoutput->parameterVals[1];
00261     (*tiles)[cnt].params.psi0 = NDoutput->coordinateVals[0];
00262     (*tiles)[cnt].params.psi3 = NDoutput->coordinateVals[1];
00263     (*tiles)[cnt].params.beta = NDoutput->coordinateVals[2];
00264     /*(*tiles)[cnt].params.eta = tmplt->eta;
00265     (*tiles)[cnt].params.chirpMass = tmplt->chirpMass;*/
00266     NDoutput = NDoutput->next;
00267   } /* for(tmplt...) */
00268                                                                                                                                                 
00269   DETATCHSTATUSPTR(status);
00270   RETURN(status);
00271 }
00272 
00273 void
00274 LALInspiralSpinBankBoundary(
00275    LALStatus            *status,
00276    NDTemplateBankInput  *input,
00277    NDTemplateBankOutput *output,
00278    INT2                 *testFlag
00279    )
00280   {
00281   REAL4 f0, mass, eta, m1, m2, m1Min, m1Max, m2Min, m2Max, betaMax, x, y, z;
00282   x = output->coordinateVals[0];
00283   y = output->coordinateVals[1];
00284   z = output->coordinateVals[2];
00285   f0 = input->f0;
00286   m1Min = input->minParameters[0];
00287   m1Max = input->maxParameters[0];
00288   m2Min = input->minParameters[1];
00289   m2Max = input->maxParameters[1];
00290   betaMax = input->maxParameters[2];
00291   *testFlag =  1;
00292  
00293   mass = -y/x / (16.0*LAL_PI*LAL_PI*f0);
00294   eta = 16.0457 * pow( -x*x/y/y/y/y/y, 0.3333333 );
00295   if (eta > 0.25 || eta < 0)
00296     *testFlag = 0;
00297   output->parameterVals[0] = m1 = 0.5*mass* (1 + sqrt(1 - 4*eta));
00298   output->parameterVals[1] = m2 = 0.5*mass* (1 - sqrt(1 - 4*eta));
00299   if (m1 > m1Max || m1 < m1Min || m2 > m2Max || m2 < m2Min)
00300     *testFlag = 0;
00301   betaMax = 3.8*LAL_PI/29.961432 * (1+0.75*m2/m1)*(m1/m2) * pow((LAL_MTSUN_SI*100.0/mass),0.6666667);
00302   if (z > betaMax)
00303     *testFlag = 0;
00304   }
00305   
00306 
00307 
00308 
00309 
00310