kat_aux.c 152 KB
Newer Older
1
/*!
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
 * \file kat_aux.c
 * \brief Auxilliary functions for Finesse
 * 
 * @section Copyright notice
 *
 *  This file is part of the interferometer simulation Finesse
 *  http://www.gwoptics.org/finesse
 *
 *  Copyright (C) 1999 onwards Andreas Freise
 *  with parts of the code written by Daniel Brown, Paul Cochrane
 *  and Gerhard Heinzel.
 *
 *  This program is free software; you can redistribute it and/or modify it under
 *  the terms of the GNU General Public License version 3 as published
 *  by the Free Software Foundation.
 *
 *  This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
 *  without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
Daniel Brown's avatar
Daniel Brown committed
20
 *  See the GNU General Public License for more details.
21
 *
Daniel Brown's avatar
Daniel Brown committed
22
 *  You should have received a copy of the GNU General Public License along with
23 24 25
 *  this library; if not, write to the Free Software Foundation, Inc., 59 Temple Place,
 *  Suite 330, Boston, MA 02111-1307 USA
 */
26 27 28 29 30
#include "kat.h"
#include "kat_inline.c"
#include "kat_io.h"
#include "kat_aux.h"
#include <errno.h>
31

32
// this include stuff is for the function getNumCores
Daniel Brown's avatar
Daniel Brown committed
33 34
#if (OS == __WIN_32__ || OS == __WIN_64__) && defined(OSWIN)

Daniel Brown's avatar
Daniel Brown committed
35
#if defined(__CYGWIN__)
36
#undef UCHAR // formulc.h defines this already and so does windows.h, so dont do it again!
37
#include <windows.h>
38 39
#include <tlhelp32.h>
#include <unistd.h>
Daniel Brown's avatar
Daniel Brown committed
40
#else
41 42 43 44
#undef UCHAR // formulc.h defines this already and so does windows.h, so dont do it again!
#include <windows.h>
#include <sys/time.h>
#include <time.h>
Daniel Brown's avatar
Daniel Brown committed
45
#endif
46
#endif
47

Daniel Brown's avatar
Daniel Brown committed
48
#if OS==MACOS
49 50 51
#include <assert.h>
#include <CoreServices/CoreServices.h>
#include <unistd.h>
Daniel Brown's avatar
Daniel Brown committed
52

53 54
#include <sys/param.h>
#include <sys/sysctl.h>
55 56 57 58 59 60
#include <sys/types.h>
#include <mach/mach.h>
#include <mach/clock.h>
#include <mach/mach_time.h>
#endif

Daniel Brown's avatar
Daniel Brown committed
61 62
#if OS == UNIX
#include <unistd.h>	
Daniel Brown's avatar
Daniel Brown committed
63
#include <time.h>
Daniel Brown's avatar
Daniel Brown committed
64
#include <sys/sysinfo.h>
Daniel Brown's avatar
Daniel Brown committed
65 66
#endif

67
#if INCLUDE_NICSLU == 1
68
#include "nicsluc.h"
69 70
#endif

71 72
//! \brief Local array to store previously computed values for the 
//! function double fac(int n) 
Daniel Brown's avatar
Daniel Brown committed
73 74
double natLogsOfN[101] = {0};
int natLogMax; //!< Maximum index of stored natural logs
75

76
extern init_variables_t init;
77 78
extern interferometer_t inter;
extern memory_t mem;
79
extern options_t options;
80
extern local_var_t vlocal;
81

82 83 84
extern time_t now;
extern timespec_t clock_start;

85 86
extern ifo_matrix_vars_t M_ifo_car;
extern ifo_matrix_vars_t M_ifo_sig;
87

88 89 90 91 92 93
extern const complex_t complex_i; //!<  sqrt(-1) or 0 + i
extern const complex_t complex_mi; //!<  sqrt(-1) or 0 - i
extern const complex_t complex_1; //!<  1 in complex space: 1 + 0i
extern const complex_t complex_m1; //!< -1 in complex space: -1 + 0i
extern const complex_t complex_0; //!<  0 in complex space: 0 + 0i

Daniel Brown's avatar
Daniel Brown committed
94
static int num_names = 0; //!< the number of names
95 96 97
char get_path_buf[LINE_LEN]; // internal buffer used by the get_path routines
                             // pointer to this is returned for reading

98 99 100
extern int *t_s; //!< Type_of_Signal list [frequency]
extern double *f_s; //!< Frequency value list [frequency]

101 102
extern int carrier_frequencies_tuned;
extern int signal_frequencies_tuned;
103

104
#ifdef OSMAC
105
mach_timebase_info_data_t    sTimebaseInfo;
106 107
#endif

108 109
char timerNames[MAX_TIMERS][TIMER_NAME_LEN];

110
timespec_t startTimes[MAX_TIMERS];
Daniel Brown's avatar
Daniel Brown committed
111
int currTimer = -1;
112
;
113 114 115 116 117

int getTime(timespec_t *tp){
#if defined(__APPLE__) && defined(__MACH__)
    kern_return_t retval = KERN_SUCCESS;
    mach_timespec_t mts;
118 119
    clock_serv_t cclock;
    
120 121 122
    host_get_clock_service(mach_host_self(), CALENDAR_CLOCK, &cclock);
    retval = clock_get_time(cclock, &mts);
    mach_port_deallocate(mach_task_self(), cclock);
123
    
124 125
    tp->tv_sec = mts.tv_sec;
    tp->tv_nsec = mts.tv_nsec;
126
    
127 128 129 130 131 132
    return retval;
#else
    return clock_gettime(CLOCK_REALTIME, tp);
#endif
}

133 134 135 136 137 138 139 140 141
int getTimeHighRes(timespec_t *tp){
#if defined(__APPLE__) && defined(__MACH__)
    tp->mach = mach_absolute_time();
    return 0;
#else
    return clock_gettime(CLOCK_REALTIME, tp);
#endif
}

142
#if defined(OSWIN) && !defined(__CYGWIN__) 
143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163
#include <string.h>

// Taken from http://www.geekdroppings.com/2014/03/16/cross-compiling-with-mingw-and-crosstool-ng/
// Added for support in mingw. This ought to be checked and enabled with autotools.
const char *strcasestr(const char *s1, const char *s2)
{
 // if either pointer is null
 if (s1 == 0 || s2 == 0)
  return 0;
 // the length of the needle
 size_t n = strlen(s2);
 // iterate through the string
 while(*s1)
 // if the compare which is case insensitive is a match, return the pointer
 if(!strnicmp(s1++,s2,n))
  return (s1-1);
 // no match was found
 return 0;
}
#endif

164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221
#ifdef __CYGWIN__

static struct pid {
	struct pid *next;
	FILE *fp;
	pid_t pid;
} *pidlist; 
	
FILE * popen2(const char *program, const char *type)
{
	struct pid *cur;
	FILE *iop;
	int pdes[2], pid;

       if ((*type != 'r' && *type != 'w')
	   || (type[1]
	       && (type[2] || (type[1] != 'b' && type[1] != 't'))
			       )) {
		errno = EINVAL;
		return (NULL);
	}

	if ((cur = malloc(sizeof(struct pid))) == NULL)
		return (NULL);

	if (pipe(pdes) < 0) {
		free(cur);
		return (NULL);
	}

	switch (pid = fork()) {
	case -1:			/* Error. */
		(void)close(pdes[0]);
		(void)close(pdes[1]);
		free(cur);
		return (NULL);
		/* NOTREACHED */
	case 0:				/* Child. */
		if (*type == 'r') {
			if (pdes[1] != STDOUT_FILENO) {
				(void)dup2(pdes[1], STDOUT_FILENO);
				(void)close(pdes[1]);
			}
			(void) close(pdes[0]);
		} else {
			if (pdes[0] != STDIN_FILENO) {
				(void)dup2(pdes[0], STDIN_FILENO);
				(void)close(pdes[0]);
			}
			(void)close(pdes[1]);
		}
        
        // search for program
		int err = execlp("cmd.exe", "cmd.exe", "/c", program, NULL);
        
        if(inter.debug)
            fprintf(stderr, "Could not find 'cmd' for popen in PATH '%s': %i err: %s \n",program,err, strerror(errno));
        
222 223
        // typically on any windows system cmd should be found
        err = execl("sh", "sh", "-c", program, NULL);
224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280
        
        if(inter.debug)
            fprintf(stderr, "Could not find 'sh' for popen '%s': %i err: %s \n",program,err, strerror(errno));
        
		err = execlp("sh", "sh", "-c", program, NULL);
        
        if(inter.debug)
            fprintf(stderr, "Could not find 'sh' for popen in PATH '%s': %i err: %s \n",program,err, strerror(errno));
        
        
        
        /// ddb - add error call instead of just exiting
        bug_error("Could not find sh.exe for popen call. Program: %s", program);
		//_exit(127);
		/* NOTREACHED */
	}

	/* Parent; assume fdopen can't fail. */
	if (*type == 'r') {
		iop = fdopen(pdes[0], type);
		(void)close(pdes[1]);
	} else {
		iop = fdopen(pdes[1], type);
		(void)close(pdes[0]);
	}

	/* Link into list of file descriptors. */
	cur->fp = iop;
	cur->pid =  pid;
	cur->next = pidlist;
	pidlist = cur;

	return (iop);
}

/*
 * pclose --
 *	Pclose returns -1 if stream is not associated with a `popened' command,
 *	if already `pclosed', or waitpid returns an error.
 */

int pclose2(FILE *iop)
{
	register struct pid *cur, *last;
	int pstat;
	pid_t pid;

	(void)fclose(iop);

	/* Find the appropriate file pointer. */
	for (last = NULL, cur = pidlist; cur; last = cur, cur = cur->next)
		if (cur->fp == iop)
			break;
	if (cur == NULL)
		return (-1);

	do {
281
            pid = waitpid(cur->pid, &pstat, 0);
282 283 284 285
	} while (pid == -1 && errno == EINTR);

	/* Remove the entry from the linked list. */
	if (last == NULL)
286
            pidlist = cur->next;
287
	else
288 289
            last->next = cur->next;
        
290 291 292 293 294 295 296 297
	free(cur);
		
	return (pid == -1 ? -1 : pstat);
}

#endif


298 299 300 301 302 303 304 305 306 307 308
/**
 * Returns a corrected path in a unix format, this is really for use on Windows
 * systems where the styling is different. If used on a unix based system then
 * the same path will be returned. This only works with Cygwin as the cygpath 
 * utility is called.
 * 
 * 
 * @param input Path to be changed
 * @param output predefined char array to put new path in
 * @return 
 */
309
void get_path_io(const char* input, char* output, bool forUnix) {
310

311
#ifdef __CYGWIN__    
312 313 314
    // silence compiler wanring
    (void) forUnix;
    
315 316 317 318
    char buf[LINE_LEN] = {0};
    char out[LINE_LEN] = {0};
    FILE *pipe;
    
319
    if(forUnix)
320
        sprintf(buf, "cygpath.exe --unix '%s' ", input);
321
    else
322
        sprintf(buf, "cygpath.exe --windows '%s' ", input);
323
    
324
    if ((pipe = popen2(buf, "rt")) == NULL)
325
        bug_error("Failed to execute popen command in get_path_io: '%s'  error: '%s'", buf,strerror(errno));
326
    
327
    if(fgets(out, LINE_LEN, pipe) == NULL)
328
        bug_error("Failed to read popen buffer in get_path_io: '%s'", buf);
329
    
330 331
    // removing trailing new line character that fgets seems to add
    size_t ln = strlen(out) - 1;
332
    
333 334
    if (out[ln] == '\n')
        out[ln] = '\0';
335
    
336 337 338 339
    pclose2(pipe);
    
    if(input != output)
        strcpy(output, out);
340
#else
341 342
    if(input != output)
        strcpy(output, input);
343
#endif
344
    (void) forUnix; //silence compiler warning.
345 346 347
}

/**
348 349 350
 * If called on anything but Cygwin compiled version then this just returns the same
 * input pointer.
 * 
351 352 353 354 355 356 357 358
 * Calls get_path_io but uses internal buffer. This should only
 * be used for instaneous reuse of the returned pointer. It is not gauranteed that
 * the value will be the same later on as get_path_i might be called again elsewhere
 * for a different path.
 * 
 * @param input path to check
 * @return converted path in unix format
 */
359
const char* get_path_i(const char* input, bool unixStyle) {
360
#ifdef __CYGWIN__
361
    get_path_io(input, get_path_buf, unixStyle);
362
    return get_path_buf;
363 364 365
#else
    return input;
#endif
366
    (void) unixStyle; // silence compiler warning
367
}
368

369 370 371 372 373 374 375 376 377 378 379 380 381
/**
 * A wrapper for System(const char* command) that when used with Cygwin will use
 * sh.exe that is present in the local directory rather than looking in path or 
 * elsewhere.
 * 
 * @param command
 * @return 
 */
int system_call(const char* command){
#ifdef __CYGWIN__
    FILE *pipe;
    
    if ((pipe = popen2(command, "rt")) == NULL)
382
        gerror("Failed to execute system_call: '%s'  error: '%s'\n", command, strerror(errno));
383 384 385 386 387 388 389 390
    
    pclose2(pipe);
    
    return errno;
#else 
    int err = system(command);
    
    if(err != 0)
391
        gerror("Failed to execute system_call: '%s'  error: '%i'\n", command, errno);
392 393 394
    
    return err;
#endif
395
 
396 397
}

398
#if INCLUDE_NICSLU == 1
399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444
int handle_NICSLU(int ret){
    switch(ret){
        case -1:
            gerror("General failure in NICSLU\n");
        case -2:
            gerror("Argument failure in NICSLU\n");
        case -3:
            gerror("Memory overflow in NICSLU\n");    
        case -4:
            gerror("Cannot open file in NICSLU\n");
        case -5:
            gerror("NICSLU matrix is structurally singular\n");
        case -6:
            warn("NICSLU matrix is numerically singular\n"); break;
        case -7:
            gerror("NICSLU matrix structure is invalid\n");
        case -8:
            gerror("NICSLU matrix entry duplicated\n");
        case -9:
            gerror("NICSLU threads not initialised\n");
        case -10:
            gerror("NICSLU matrix not created\n");
        case -11:
            gerror("NICSLU thread scheduler not intialised\n");
        case -12:
            gerror("NICSLU cannot create single thread\n");
        case -13:
            gerror("NICSLU thread cannot be created\n");
        case -14:
            gerror("NICSLU matrix not analysed\n");
        case -15:
            gerror("NICSLU matrix not factorised\n");
        case 1:
            return NICSLU_USE_SEQUENTIAL_FACTORIZATION;
        case 2:
            gerror("NICSLU threads cannot be pinned to core\n");
            
        case 0:
        default:
            return NICS_OK;
    }
    
    return NICS_OK;
}
#endif

Daniel Brown's avatar
Daniel Brown committed
445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482
/**
 * For a given frequency the carrier field frequency list is searched and the
 * equivalent frequency object is returned. Returns NULL if not found.
 * 
 * @param frequency
 * @return 
 */
frequency_t* get_carrier_frequency(double frequency){
    int f;
    
    for(f=0; f<inter.num_carrier_freqs; f++){
        if(eq(frequency, inter.carrier_f_list[f]->f)){
            return inter.carrier_f_list[f];
        }
    }
    
    return NULL;
}

/**
 * For a given frequency the signal field frequency list is searched and the
 * equivalent frequency object is returned. Returns NULL if not found.
 * 
 * @param frequency
 * @return 
 */
frequency_t* get_signal_frequency(double frequency){
    int f;
    
    for(f=0; f<inter.num_signal_freqs; f++){
        if(eq(frequency, inter.signal_f_list[f]->f)){
            return inter.signal_f_list[f];
        }
    }
    
    return NULL;
}

483
//! The type of signal of the fsig command
Daniel Brown's avatar
Daniel Brown committed
484

485
/*!
486
 * \param s the input string to determine the signal type
487
 *
488
 * \see Test_get_signal_type()
489
 */
Daniel Brown's avatar
Daniel Brown committed
490
int get_signal_type(char *signal_type_string) {
491 492
    assert(signal_type_string != NULL);
    
493
    if (string_matches_exactly(signal_type_string, "phase")) {
Daniel Brown's avatar
Daniel Brown committed
494
        return SIG_PHS;
495 496
    } else if (string_matches_exactly(signal_type_string, "z")) {
        return SIG_Z;
497 498
    } else if (string_matches_exactly(signal_type_string, "phi")) {
        return SIG_PHS;
499
    } else if (string_matches_exactly(signal_type_string, "phs")) {
500
        return SIG_PHS;
501
    } else if (string_matches_exactly(signal_type_string, "amp")) {
Daniel Brown's avatar
Daniel Brown committed
502
        return SIG_AMP;
503
    } else if (string_matches_exactly(signal_type_string, "freq")) {
Daniel Brown's avatar
Daniel Brown committed
504
        return SIG_FRQ;
505
    } else if (string_matches_exactly(signal_type_string, "frq")) {
Daniel Brown's avatar
Daniel Brown committed
506
        return SIG_FRQ;
507
    } else if (string_matches_exactly(signal_type_string, "xbeta") ||
508 509
            string_matches_exactly(signal_type_string, "rx") ||
            string_matches_exactly(signal_type_string, "yaw")) {
Daniel Brown's avatar
Daniel Brown committed
510
        return SIG_X;
511
    } else if (string_matches_exactly(signal_type_string, "ybeta") ||
512 513
            string_matches_exactly(signal_type_string, "ry") ||
            string_matches_exactly(signal_type_string, "pitch")) {
Daniel Brown's avatar
Daniel Brown committed
514
        return SIG_Y;
515
    } else if (string_matches_exactly(signal_type_string, "p")) {
Daniel Brown's avatar
Daniel Brown committed
516
        return SIG_PHS;
517
    } else if (string_matches_exactly(signal_type_string, "a")) {
Daniel Brown's avatar
Daniel Brown committed
518
        return SIG_AMP;
519
    } else if (string_matches_exactly(signal_type_string, "f")) {
Daniel Brown's avatar
Daniel Brown committed
520
        return SIG_FRQ;
521
    } else if (string_matches_exactly(signal_type_string, "x")) {
Daniel Brown's avatar
Daniel Brown committed
522
        return SIG_X;
523
    } else if (string_matches_exactly(signal_type_string, "y")) {
Daniel Brown's avatar
Daniel Brown committed
524
        return SIG_Y;
525 526
    } else if (string_matches_exactly(signal_type_string, "Rc")) {
        return SIG_RC;
Daniel Brown's avatar
Daniel Brown committed
527 528
    } else if (string_matches_exactly(signal_type_string, "sagnac")) {
        return SIG_SAGNAC;
529
    } else if (signal_type_string[0] == 's' && isdigit(signal_type_string[1])) {
530 531
        // check if we are trying to drive a surface mode
        return SIG_SURF;
532 533 534 535 536 537
    } else if (string_matches_exactly(signal_type_string, "Fz")) {
        return SIG_FZ;
    } else if (string_matches_exactly(signal_type_string, "Frx") || string_matches_exactly(signal_type_string, "Fyaw")) {
        return SIG_FRX;
    } else if (string_matches_exactly(signal_type_string, "Fry") || string_matches_exactly(signal_type_string, "Fpitch")) {
        return SIG_FRY;
538
    } else if ((signal_type_string[0] == 'F' && signal_type_string[1] == 's') && isdigit(signal_type_string[2])) {
539 540
        // check if we are trying to drive a surface mode
        return SIG_FSURF;
Daniel Brown's avatar
Daniel Brown committed
541 542
    }
    return (NOT_FOUND);
543 544
}

545 546
char indent[MAX_TIMERS] = "";
FILE *timeFile = NULL;
Daniel Brown's avatar
Daniel Brown committed
547
bool enableTiming = false;
548

549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584
int motion_string_to_index(const char *v, double M, double Ix, double Iy, int num_surf_motion){
    
    if(strcasecmp("z", v) == 0){
        if (M == 0.0) {
            return NOT_FOUND;
        } else
            return 0;
    } else if(strcasecmp("rx", v) == 0 || strcasecmp("yaw", v) == 0){
        if(Ix == 0.0) {
            return NOT_FOUND;
        } else
            return (M==0) ? 0 : 1;
    } else if(strcasecmp("ry", v) == 0 || strcasecmp("pitch", v) == 0){
        if(Iy == 0.0) {
            return NOT_FOUND;
        } else if (M==0 && Ix == 0) return 0;
        else if ((M != 0) ^ (Ix != 0)) return 1;
        else if ((M != 0) && (Ix != 0)) return 2;
        else return NOT_FOUND;
    } else if(strlen(v) >= 2 && ('s' == v[0] || 'S' == v[0])){
        int i=0;
        
        if(M!=0) i++;
        if(Ix!=0) i++;
        if(Iy!=0) i++;
        
        int ix = atoi(&v[1]);
        
        if(ix < num_surf_motion)
            return i + ix;
        else
            return NOT_FOUND;
    }
    
    return NOT_FOUND;
}
585

586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660
int motion_string_to_mirror_tf(const char *v, mirror_t *m, transfer_func_t **tf){
    double M = m->mass;
    double Ix = m->Ix;
    double Iy = m->Iy;
    *tf = NULL;
    
    if(strcasecmp("z", v) == 0){
        *tf = m->long_tf;
        if (M == 0.0) {
            return NOT_FOUND;
        } else
            return 0;
    } else if(strcasecmp("rx", v) == 0 || strcasecmp("yaw", v) == 0){
        *tf = m->rot_x_tf;
        if(Ix == 0.0) {
            return NOT_FOUND;
        } else
            return (M==0) ? 0 : 1;
    } else if(strcasecmp("ry", v) == 0 || strcasecmp("pitch", v) == 0){
        *tf = m->rot_y_tf;
        if(Iy == 0.0) {
            return NOT_FOUND;
        } else if (M==0 && Ix == 0) return 0;
        else if ((M != 0) ^ (Ix != 0)) return 1;
        else if ((M != 0) && (Ix != 0)) return 2;
        else return NOT_FOUND;
    } else if(strlen(v) >= 2 && ('s' == v[0] || 'S' == v[0])){
        int i=0;
        
        if(M!=0) i++;
        if(Ix!=0) i++;
        if(Iy!=0) i++;
        
        int ix = atoi(&v[1]);
        
        if(ix < m->num_surface_motions){
            *tf = m->surface_motions_tf[ix];
            return i + ix;
        } else
            return NOT_FOUND;
    }
    
    return NOT_FOUND;
}

int motion_string_to_bs_tf(const char *v, beamsplitter_t *m, transfer_func_t **tf){
    double M = m->mass;
    double Ix = m->Ix;
    double Iy = m->Iy;
    
    if(strcasecmp("z", v) == 0){
        *tf = m->long_tf;
        if (M == 0.0) {
            return NOT_FOUND;
        } else
            return 0;
    } else if(strcasecmp("rx", v) == 0 || strcasecmp("yaw", v) == 0){
        *tf = m->rot_x_tf;
        if(Ix == 0.0) {
            return NOT_FOUND;
        } else
            return (M==0) ? 0 : 1;
    } else if(strcasecmp("ry", v) == 0 || strcasecmp("pitch", v) == 0){
        *tf = m->rot_y_tf;
        if(Iy == 0.0) {
            return NOT_FOUND;
        } else if (M==0 && Ix == 0) return 0;
        else if ((M != 0) ^ (Ix != 0)) return 1;
        else if ((M != 0) && (Ix != 0)) return 2;
        else return NOT_FOUND;
    }
    
    return NOT_FOUND;
}

661 662 663 664 665 666 667 668
/**
 * Returns index of motion for a given string input
 */
int parse_motion_value(motion_out_t *mout, double M, double Ix, double Iy, int num_surf_motion){
    const char *v = mout->motion_name;
    const char *cname = get_comp_name2(mout->component_type, mout->component_index);
    
    if(strcasecmp("z", v) == 0){
669 670 671 672 673 674 675
        if (M == 0.0) {
            if(!(mout->warnings & LONG_WARN)) {
                warn("The longitudinal motion does not exist for %s (detector %s)\n", cname, mout->name);
                mout->warnings |= LONG_WARN;
            }
            return NOT_FOUND;
        } else
676 677
            return 0;
    } else if(strcasecmp("rx", v) == 0 || strcasecmp("yaw", v) == 0){
678 679 680 681 682 683 684 685
        if(Ix == 0.0) {
            if(!(mout->warnings & ROTX_WARN)) {
                warn("The yaw motion does not exist for %s (detector %s)\n", cname, mout->name);
                mout->warnings |= ROTX_WARN;
            }
            
            return NOT_FOUND;
        } else
686 687
            return (M==0) ? 0 : 1;
    } else if(strcasecmp("ry", v) == 0 || strcasecmp("pitch", v) == 0){
688 689 690 691 692 693 694 695 696
        if(Iy == 0.0) {
            if(!(mout->warnings & ROTY_WARN)) {
                warn("The pitch motion does not exist for %s (detector %s)\n", cname, mout->name);
                mout->warnings |= ROTY_WARN;
            }
            
            warn("The pitch motion does not exist for %s (detector %s)\n", cname, mout->name);
            return NOT_FOUND;
        } else if (M==0 && Ix == 0) return 0;
Daniel Brown's avatar
Daniel Brown committed
697 698
        else if ((M != 0) ^ (Ix != 0)) return 1;
        else if ((M != 0) && (Ix != 0)) return 2;
699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717
        else return NOT_FOUND;
    } else if(strlen(v) >= 2 && ('s' == v[0] || 'S' == v[0])){
        int i=0;
        
        if(M!=0) i++;
        if(Ix!=0) i++;
        if(Iy!=0) i++;
        
        int ix = atoi(&v[1]);
        
        if(ix >= num_surf_motion) gerror("The surface motion '%s' does not exist at component %s for detector %s\n", v, cname, mout->name);
        
        return i + ix;
    } else
        gerror("The motion '%s' does not exist for component %s (detector %s)\n", v, cname, mout->name);
    
    return NOT_FOUND;
}

718

719 720 721 722 723 724 725 726 727
int is_minimizer_tuned(void *var, void **rtn, int *type){
    
    if(inter.minimizer == NULL)
        return 0;
    
    if(inter.minimizer->variable != var)
        return 0;
    
    
Daniel Brown's avatar
Daniel Brown committed
728 729
    if(rtn)*rtn = &inter.minimizer;
    if(type)*type = MINIMIZER;
730 731 732 733
          
    return 1;
}

734 735 736 737 738 739 740 741 742
/**
 * Determines if a double variable is being tuned by any of the xaxes. If so it
 * returns which one.
 * 
 * @param var
 * @param rtn
 * @param type
 * @return 
 */
743 744 745
int is_xaxis_tuned(void *var, void **rtn, int *type){
    
    if(var == inter.x1.xaxis){
Daniel Brown's avatar
Daniel Brown committed
746 747
        if(rtn)*rtn = &inter.x1;
        if(type)*type = X1;
748 749
        return 1;
    } else if(var == inter.x2.xaxis){
Daniel Brown's avatar
Daniel Brown committed
750 751
        if(rtn)*rtn = &inter.x2;
        if(type)*type = X2;
752 753
        return 1;
    } else if(var == inter.x3.xaxis){
Daniel Brown's avatar
Daniel Brown committed
754 755
        if(rtn)*rtn = &inter.x3;
        if(type)*type = X3;
756 757 758
        return 1;
    }
    
Daniel Brown's avatar
Daniel Brown committed
759 760
    if(type)*type = -1;
    if(rtn)*rtn = NULL;
761 762 763
    return 0;
}

764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807
void bs_get_nr(beamsplitter_t *bs, double *nr1, double *nr2){
    assert(bs  != NULL);
    assert(nr1 != NULL);
    assert(nr2 != NULL);
    
    node_t *node1 = &inter.node_list[bs->node1_index];
    node_t *node2 = &inter.node_list[bs->node2_index];
    node_t *node3 = &inter.node_list[bs->node3_index];
    node_t *node4 = &inter.node_list[bs->node4_index];
    
    assert(node1 != NULL);
    assert(node2 != NULL);
    assert(node3 != NULL);
    assert(node4 != NULL);  
    
    if(!node1->gnd_node && !node2->gnd_node){
        if((*node1->n != *node2->n)){
            gerror("Refractive indices at node 1 and node 2 do not match for %s\n", bs->name);
        }
    }
    
    if(!node3->gnd_node && !node4->gnd_node){
        if((*node3->n != *node4->n)){
            gerror("Refractive indices at node 3 and node 4 do not match for %s\n", bs->name);
        } 
    }
    
    if(!node1->gnd_node){
        *nr1 = *node1->n;
    } else if(!node2->gnd_node){
        *nr1 = *node2->n;
    } else {
        *nr1 = 1;
    }
    
    if(!node3->gnd_node) {
        *nr2 = *node3->n;
    } else if(!node4->gnd_node) {
        *nr2 = *node4->n;
    } else {
        *nr2 = 1;
    }
}

808 809 810 811 812
int is_put_tuned(void *var, void **put, int *type) {
    int i;
    
    for(i=0; i < inter.num_put_cmds; i++){
        if(inter.put_list[i].target == var){
813 814
            if(put)  *put = &inter.put_list[i];
            if(type) *type = PUT;
815 816 817 818
            return 1;
        }
    }
    
819 820 821
    if(type) *type = -1;
    if(put)  *put = NULL;
    
822 823 824
    return 0; 
}

825 826 827 828 829 830 831 832
/**
 * Checks to see if a variable is being tuned by either a:
 *  - minimizer
 *  - xaxis, x2axis
 *  - put
 * 
 *  source: pointer to double value that the tuning value comes from
 *  rtn   : is set to a pointer to the object that is tuning.
Daniel Brown's avatar
Daniel Brown committed
833
 *  type  : is the type of thing doing the tuning
834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849
 * 
 * @return True or False
 */
bool is_tuned(double *var, double **source, void **tuner, int *type){
    bool rtn = false;
    
    if(!rtn){
        if(is_minimizer_tuned(var, tuner, type)){
            assert(true == false && "Not implemented yet");
            rtn = true;
        }
    }
    
    if(!rtn){
        if(is_xaxis_tuned(var, tuner, type)) {
            axis_t *axis = *tuner;
Daniel Brown's avatar
Daniel Brown committed
850
            if(source) *source = axis->xaxis;
851 852 853 854 855 856 857
            rtn = true;
        }
    }
    
    if(!rtn){
        if(is_put_tuned(var, tuner, type)){
            put_command_t *put = *tuner;
Daniel Brown's avatar
Daniel Brown committed
858
            if(source) *source = put->source;
859 860 861 862 863 864 865 866 867
            rtn = true;
        }
    }
    
    
    
    return rtn;
}

868 869 870 871 872 873 874 875 876 877 878 879
timespec_t diff(timespec_t start, timespec_t end){
	timespec_t temp;
    
	if ((end.tv_nsec-start.tv_nsec)<0) {
		temp.tv_sec = end.tv_sec-start.tv_sec-1;
		temp.tv_nsec = 1000000000+end.tv_nsec-start.tv_nsec;
	} else {
		temp.tv_sec = end.tv_sec-start.tv_sec;
		temp.tv_nsec = end.tv_nsec-start.tv_nsec;
	}
    
	return temp;
880
}
881

Daniel Brown's avatar
Daniel Brown committed
882
int startTimer(const char* name){
Daniel Brown's avatar
Daniel Brown committed
883
    if(!enableTiming) return -1;
884
    
Daniel Brown's avatar
Daniel Brown committed
885 886 887 888 889
    if(timeFile == NULL){
        char buf[LINE_LEN];
        sprintf(buf, "%s.perf", inter.basename);
        timeFile = fopen(buf,"w");
    }
890 891 892 893 894 895
    
    currTimer++;
    
    if(currTimer >= MAX_TIMERS)
        bug_error("Tried to create more timers than was allowed: %i out of %i %s",currTimer, MAX_TIMERS, name);
    
896 897
    getTimeHighRes(&startTimes[currTimer]);
    
898 899 900 901 902 903
    strcpy(timerNames[currTimer], name); // record what we are timing
    strcat(indent, " "); // increase indent for visual separation of nested timers
    
    return currTimer;
}

Daniel Brown's avatar
Daniel Brown committed
904
void endTimer(int id){
905 906 907 908 909 910 911
    if(!enableTiming) return;
    
    if(currTimer >= MAX_TIMERS)
        bug_error("number of timers was greater than MAX_TIMERS on ending timer %i", id);
    
    if(id != currTimer)
        bug_error("Tried to stop timer %i (%s) before timer %i (%s)", id, timerNames[id], currTimer, timerNames[currTimer]);
912
    
913 914
    timespec_t endt;
    getTimeHighRes(&endt);
915
    
916 917 918 919 920 921 922 923
#if defined(__APPLE__) && defined(__MACH__)
    uint64_t elapsedNano = (endt.mach - startTimes[currTimer].mach) * sTimebaseInfo.numer / sTimebaseInfo.denom;
    //printf("%llu\n", elapsedNano);
    fprintf(timeFile,"%s%s %lu %llu\n", indent, timerNames[currTimer], 0L, elapsedNano);
#else
    timespec_t ets = diff(startTimes[currTimer], endt);
    fprintf(timeFile,"%s%s %lu %lu\n", indent, timerNames[currTimer], ets.tv_sec, ets.tv_nsec);
#endif
924 925 926 927 928
    
    indent[currTimer] = '\0';
    
    currTimer--;
}
929

930 931 932 933 934 935 936 937 938 939

void closeTimerFile(){
    if(!enableTiming) return;
    
    if(timeFile != NULL){
        fflush(timeFile);
        fclose(timeFile);
    }
}

940 941 942
int _find_frequency_in_list(frequency_t *f, frequency_node_t **found){
    assert(found != NULL);
    assert(f != NULL);
943
    
944
    frequency_node_t *curr = inter.frequency_linked_list.head;
945
    
946
    while(curr != NULL){
947
        if(curr->me == f){
948
            *found = curr;
949 950
            return 1;
        }
951 952
        
        curr = curr->next;
953 954
    }
    
955
    *found = NULL;
956 957 958
    return 0;
}

959 960 961 962 963 964 965 966 967 968 969 970 971 972 973 974 975
/**
 * Given a pointer to a double, this function will check to see if this variable
 * is being put or tuned by an xaxis command.
 * 
 * @param var pointer to variable to check
 * @return true if being tuned, false if not
 */
bool is_var_tuned(double *var) {
    if(!inter.noxaxis){
        if(var == inter.x1.xaxis) return true;
        if(var == inter.x2.xaxis) return true;
        if(inter.x3_axis_is_set && var == inter.x3.xaxis) return true;
    }
    
    int i;
    
    for(i=0; i<inter.num_put_cmds; i++){
976
        if(inter.put_list[i].target == var) return true;
977 978
    }
    
979 980 981
    if(inter.minimizer != NULL && inter.minimizer->variable == var)
        return true;
    
982 983 984
    return false;
}

985 986
void _delete_frequency_from_list(frequency_node_t *node){
    assert(node != NULL);
987 988
    node->me->removed = true;
    node->me->node = NULL;
989 990 991 992 993 994 995 996 997 998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010 1011 1012 1013 1014 1015 1016
    
    if(node->prev != NULL && node->next != NULL){
        // if there are nodes before and after this...
        node->prev->next = node->next;
        node->next->prev = node->prev;
    } else if(node->prev == NULL && node->next != NULL){
        // if there is no node before, but one after..
        node->next->prev = NULL;
    } else if(node->prev != NULL && node->next == NULL){
        // if there is no node after, but one before..
        node->prev->next = NULL;
    }
    
    if(inter.frequency_linked_list.head == node)
        inter.frequency_linked_list.head = node->next;
    
    if(inter.frequency_linked_list.end == node)
        inter.frequency_linked_list.end = node->prev;
    
    free(node);
    
    inter.frequency_linked_list.count--;
}

void _insert_frequency_to_list(frequency_node_t *insert_after, frequency_t *freq){
    assert(freq != NULL);
    
    frequency_node_t *node = (frequency_node_t*) calloc(1, sizeof(frequency_node_t));
1017
    node->me = freq;
1018
    freq->node = node;
1019 1020 1021 1022
    
    assert(node != NULL);
    assert(node->next == NULL);
    assert(node->prev == NULL);
1023
    assert(node->me != NULL);
1024
    
1025 1026 1027 1028 1029 1030 1031 1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049 1050 1051 1052
    frequency_linked_list_t *list = &inter.frequency_linked_list;
    
    if(list->head == NULL){
        // we have an empty list
        assert(list->count == 0);
        
        list->head = node;
        list->end = node;
    } else {
        
        // if we haven't specified where to add to the list put it at the end
        if(insert_after == NULL)
            insert_after = list->end;
        
        // if there is another item next
        if(insert_after->next != NULL){
            insert_after->next->prev = node;
            node->next = insert_after->next;
        }
        
        insert_after->next = node;
        node->prev = insert_after;
                
        // if we are inserting at the end update the end node
        if(insert_after == list->end)
            list->end = node;
        
    }
1053
    
1054 1055 1056
    list->count++;
}

1057
void print_used_frequency_list(){
1058
    message("** Total Frequency List (Difference from base Frequency %E Hz):\n", inter.f0);
1059 1060 1061 1062
    
    int i, j;
    
    for(i=0; i<inter.num_carrier_freqs; i++){
1063
        frequency_t *fc = inter.carrier_f_list[i];
1064
        
Daniel Brown's avatar
Daniel Brown committed
1065
        message("   %i  %s  %16.16gHz  tuned=%i  ", fc->index, fc->name, fc->f, fc->isTuned);
1066

1067
        switch(fc->type){
1068
            case CARRIER_FIELD:
1069
                message("Carrier\n");
1070 1071
                break;
            case MODULATOR_FIELD:
1072
                message("Modulation\n");
1073 1074
                break;
            case SIGNAL_FIELD:
1075
                message("Signal\n");
1076
                break;
1077
            case USER_FIELD:
1078
                message("User\n");
1079
                break;
1080
        }
1081 1082
        
        for(j=0; j<inter.num_signal_freqs; j++){
1083
            frequency_t *fs = inter.signal_f_list[j];
1084 1085
            
            if(fs->carrier == fc)
Daniel Brown's avatar
Daniel Brown committed
1086
                message("   %i  %s  %16.16gHz  tuned=%i  Signal\n", fs->index, fs->name, fs->f, fs->isTuned);
1087
        }
1088 1089 1090
    }
}

1091 1092
void reduce_frequency_linked_list(){
    frequency_node_t *curr = inter.frequency_linked_list.head, *other = NULL;
1093 1094
    int i;
    
1095 1096 1097 1098 1099 1100 1101
    while(curr != NULL){
        // if there are no other frequencies stop as there is nothing
        // more we can remove anyway
        if(curr->next == NULL)
            break;

        // if this frequency is tuned then its value will change and we
1102
        // can't reduce it...
1103
        if(!curr->me->isTuned){
1104
            
1105
            other = inter.frequency_linked_list.head;
1106

1107 1108
            // loop through all other frequencies we haven't looked through 
            // yet and see if any of them can be removed
1109
            while(other != NULL) {
1110
                int delete_other = 0;
1111

1112
                if(!(other->me->isTuned || other == curr)){
1113 1114

                    // check if this other frequency has the same value
1115
                    if(other->me->f == curr->me->f){
1116 1117 1118
                        // as carrier fields will always be larger than signal or modulator
                        // sidebands keep them over the others. Then modulator over signals

1119
                        switch(curr->me->type){
1120 1121 1122 1123
                            case CARRIER_FIELD:
                                // if the current frequency is a carrier then delete whatever
                                // duplicate comes later
                                delete_other = 1;
1124 1125 1126
                                
                                // if we are