/*! * \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. * See the GNU General Public License for more details. * * You should have received a copy of the GNU General Public License along with * this library; if not, write to the Free Software Foundation, Inc., 59 Temple Place, * Suite 330, Boston, MA 02111-1307 USA */ #include "kat.h" #include "kat_inline.c" #include "kat_io.h" #include "kat_aux.h" #include // this include stuff is for the function getNumCores #if (OS == __WIN_32__ || OS == __WIN_64__) && defined(OSWIN) #if defined(__CYGWIN__) #undef UCHAR // formulc.h defines this already and so does windows.h, so dont do it again! #include #include #include #else #undef UCHAR // formulc.h defines this already and so does windows.h, so dont do it again! #include #include #include #endif #endif #if OS==MACOS #include #include #include #include #include #include #include #include #include #endif #if OS == UNIX #include #include #include #endif #if INCLUDE_NICSLU == 1 #include "nicsluc.h" #endif //! \brief Local array to store previously computed values for the //! function double fac(int n) double natLogsOfN[101] = {0}; int natLogMax; //!< Maximum index of stored natural logs extern init_variables_t init; extern interferometer_t inter; extern memory_t mem; extern options_t options; extern local_var_t vlocal; extern time_t now; extern timespec_t clock_start; extern ifo_matrix_vars_t M_ifo_car; extern ifo_matrix_vars_t M_ifo_sig; 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 static int num_names = 0; //!< the number of names char get_path_buf[LINE_LEN]; // internal buffer used by the get_path routines // pointer to this is returned for reading extern int *t_s; //!< Type_of_Signal list [frequency] extern double *f_s; //!< Frequency value list [frequency] extern int carrier_frequencies_tuned; extern int signal_frequencies_tuned; #ifdef OSMAC mach_timebase_info_data_t sTimebaseInfo; #endif char timerNames[MAX_TIMERS][TIMER_NAME_LEN]; timespec_t startTimes[MAX_TIMERS]; int currTimer = -1; ; int getTime(timespec_t *tp){ #if defined(__APPLE__) && defined(__MACH__) kern_return_t retval = KERN_SUCCESS; mach_timespec_t mts; clock_serv_t cclock; host_get_clock_service(mach_host_self(), CALENDAR_CLOCK, &cclock); retval = clock_get_time(cclock, &mts); mach_port_deallocate(mach_task_self(), cclock); tp->tv_sec = mts.tv_sec; tp->tv_nsec = mts.tv_nsec; return retval; #else return clock_gettime(CLOCK_REALTIME, tp); #endif } 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 } #if defined(OSWIN) && !defined(__CYGWIN__) #include // 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 #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)); // typically on any windows system cmd should be found err = execl("sh", "sh", "-c", program, NULL); 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 { pid = waitpid(cur->pid, &pstat, 0); } while (pid == -1 && errno == EINTR); /* Remove the entry from the linked list. */ if (last == NULL) pidlist = cur->next; else last->next = cur->next; free(cur); return (pid == -1 ? -1 : pstat); } #endif /** * 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 */ void get_path_io(const char* input, char* output, bool forUnix) { #ifdef __CYGWIN__ // silence compiler wanring (void) forUnix; char buf[LINE_LEN] = {0}; char out[LINE_LEN] = {0}; FILE *pipe; if(forUnix) sprintf(buf, "cygpath.exe --unix '%s' ", input); else sprintf(buf, "cygpath.exe --windows '%s' ", input); if ((pipe = popen2(buf, "rt")) == NULL) bug_error("Failed to execute popen command in get_path_io: '%s' error: '%s'", buf,strerror(errno)); if(fgets(out, LINE_LEN, pipe) == NULL) bug_error("Failed to read popen buffer in get_path_io: '%s'", buf); // removing trailing new line character that fgets seems to add size_t ln = strlen(out) - 1; if (out[ln] == '\n') out[ln] = '\0'; pclose2(pipe); if(input != output) strcpy(output, out); #else if(input != output) strcpy(output, input); #endif (void) forUnix; //silence compiler warning. } /** * If called on anything but Cygwin compiled version then this just returns the same * input pointer. * * 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 */ const char* get_path_i(const char* input, bool unixStyle) { #ifdef __CYGWIN__ get_path_io(input, get_path_buf, unixStyle); return get_path_buf; #else return input; #endif (void) unixStyle; // silence compiler warning } /** * 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) gerror("Failed to execute system_call: '%s' error: '%s'\n", command, strerror(errno)); pclose2(pipe); return errno; #else int err = system(command); if(err != 0) gerror("Failed to execute system_call: '%s' error: '%i'\n", command, errno); return err; #endif } #if INCLUDE_NICSLU == 1 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 /** * 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; ff)){ 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; ff)){ return inter.signal_f_list[f]; } } return NULL; } //! The type of signal of the fsig command /*! * \param s the input string to determine the signal type * * \see Test_get_signal_type() */ int get_signal_type(char *signal_type_string) { assert(signal_type_string != NULL); if (string_matches_exactly(signal_type_string, "phase")) { return SIG_PHS; } else if (string_matches_exactly(signal_type_string, "z")) { return SIG_Z; } else if (string_matches_exactly(signal_type_string, "phi")) { return SIG_PHS; } else if (string_matches_exactly(signal_type_string, "phs")) { return SIG_PHS; } else if (string_matches_exactly(signal_type_string, "amp")) { return SIG_AMP; } else if (string_matches_exactly(signal_type_string, "freq")) { return SIG_FRQ; } else if (string_matches_exactly(signal_type_string, "frq")) { return SIG_FRQ; } else if (string_matches_exactly(signal_type_string, "xbeta") || string_matches_exactly(signal_type_string, "rx") || string_matches_exactly(signal_type_string, "yaw")) { return SIG_X; } else if (string_matches_exactly(signal_type_string, "ybeta") || string_matches_exactly(signal_type_string, "ry") || string_matches_exactly(signal_type_string, "pitch")) { return SIG_Y; } else if (string_matches_exactly(signal_type_string, "p")) { return SIG_PHS; } else if (string_matches_exactly(signal_type_string, "a")) { return SIG_AMP; } else if (string_matches_exactly(signal_type_string, "f")) { return SIG_FRQ; } else if (string_matches_exactly(signal_type_string, "x")) { return SIG_X; } else if (string_matches_exactly(signal_type_string, "y")) { return SIG_Y; } else if (string_matches_exactly(signal_type_string, "Rc")) { return SIG_RC; } else if (string_matches_exactly(signal_type_string, "sagnac")) { return SIG_SAGNAC; } else if (signal_type_string[0] == 's' && isdigit(signal_type_string[1])) { // check if we are trying to drive a surface mode return SIG_SURF; } 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; } else if ((signal_type_string[0] == 'F' && signal_type_string[1] == 's') && isdigit(signal_type_string[2])) { // check if we are trying to drive a surface mode return SIG_FSURF; } return (NOT_FOUND); } char indent[MAX_TIMERS] = ""; FILE *timeFile = NULL; bool enableTiming = false; 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; } 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; } /** * 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){ 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 return 0; } else if(strcasecmp("rx", v) == 0 || strcasecmp("yaw", v) == 0){ 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 return (M==0) ? 0 : 1; } else if(strcasecmp("ry", v) == 0 || strcasecmp("pitch", v) == 0){ 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; 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) 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; } int is_minimizer_tuned(void *var, void **rtn, int *type){ if(inter.minimizer == NULL) return 0; if(inter.minimizer->variable != var) return 0; if(rtn)*rtn = &inter.minimizer; if(type)*type = MINIMIZER; return 1; } /** * 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 */ int is_xaxis_tuned(void *var, void **rtn, int *type){ if(var == inter.x1.xaxis){ if(rtn)*rtn = &inter.x1; if(type)*type = X1; return 1; } else if(var == inter.x2.xaxis){ if(rtn)*rtn = &inter.x2; if(type)*type = X2; return 1; } else if(var == inter.x3.xaxis){ if(rtn)*rtn = &inter.x3; if(type)*type = X3; return 1; } if(type)*type = -1; if(rtn)*rtn = NULL; return 0; } 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; } } 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){ if(put) *put = &inter.put_list[i]; if(type) *type = PUT; return 1; } } if(type) *type = -1; if(put) *put = NULL; return 0; } /** * 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. * type : is the type of thing doing the tuning * * @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; if(source) *source = axis->xaxis; rtn = true; } } if(!rtn){ if(is_put_tuned(var, tuner, type)){ put_command_t *put = *tuner; if(source) *source = put->source; rtn = true; } } return rtn; } 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; } int startTimer(const char* name){ if(!enableTiming) return -1; if(timeFile == NULL){ char buf[LINE_LEN]; sprintf(buf, "%s.perf", inter.basename); timeFile = fopen(buf,"w"); } currTimer++; if(currTimer >= MAX_TIMERS) bug_error("Tried to create more timers than was allowed: %i out of %i %s",currTimer, MAX_TIMERS, name); getTimeHighRes(&startTimes[currTimer]); strcpy(timerNames[currTimer], name); // record what we are timing strcat(indent, " "); // increase indent for visual separation of nested timers return currTimer; } void endTimer(int id){ 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]); timespec_t endt; getTimeHighRes(&endt); #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 indent[currTimer] = '\0'; currTimer--; } void closeTimerFile(){ if(!enableTiming) return; if(timeFile != NULL){ fflush(timeFile); fclose(timeFile); } } int _find_frequency_in_list(frequency_t *f, frequency_node_t **found){ assert(found != NULL); assert(f != NULL); frequency_node_t *curr = inter.frequency_linked_list.head; while(curr != NULL){ if(curr->me == f){ *found = curr; return 1; } curr = curr->next; } *found = NULL; return 0; } /** * 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; ivariable == var) return true; return false; } void _delete_frequency_from_list(frequency_node_t *node){ assert(node != NULL); node->me->removed = true; node->me->node = NULL; 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)); node->me = freq; freq->node = node; assert(node != NULL); assert(node->next == NULL); assert(node->prev == NULL); assert(node->me != NULL); 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; } list->count++; } void print_used_frequency_list(){ message("** Total Frequency List (Difference from base Frequency %E Hz):\n", inter.f0); int i, j; for(i=0; iindex, fc->name, fc->f, fc->isTuned); switch(fc->type){ case CARRIER_FIELD: message("Carrier\n"); break; case MODULATOR_FIELD: message("Modulation\n"); break; case SIGNAL_FIELD: message("Signal\n"); break; case USER_FIELD: message("User\n"); break; } for(j=0; jcarrier == fc) message(" %i %s %16.16gHz tuned=%i Signal\n", fs->index, fs->name, fs->f, fs->isTuned); } } } void reduce_frequency_linked_list(){ frequency_node_t *curr = inter.frequency_linked_list.head, *other = NULL; int i; 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 // can't reduce it... if(!curr->me->isTuned){ other = inter.frequency_linked_list.head; // loop through all other frequencies we haven't looked through // yet and see if any of them can be removed while(other != NULL) { int delete_other = 0; if(!(other->me->isTuned || other == curr)){ // check if this other frequency has the same value if(other->me->f == curr->me->f){ // as carrier fields will always be larger than signal or modulator // sidebands keep them over the others. Then modulator over signals switch(curr->me->type){ case CARRIER_FIELD: // if the current frequency is a carrier then delete whatever // duplicate comes later delete_other = 1; // if we are deleting another carrier we have to // update the frequency pointer for the light input if(other->me->type == CARRIER_FIELD){ for(i=0; ime) inter.light_in_list[i].f = curr->me; if(inter.light_in_list[i].f_entangled_zero == other->me) inter.light_in_list[i].f_entangled_zero = curr->me; } } break; case MODULATOR_FIELD: if(other->me->type == MODULATOR_FIELD || other->me->type == SIGNAL_FIELD) delete_other = 1; break; case SIGNAL_FIELD: if(other->me->type == SIGNAL_FIELD) delete_other = 1; break; } } } if(delete_other){ // move to next field then remove the older one frequency_node_t *_other = other; other = other->next; if(inter.fsig > 0 && (_other->me->type == CARRIER_FIELD || _other->me->type == MODULATOR_FIELD)) { assert(!((_other->me->sig_upper != NULL) ^ (_other->me->sig_lower != NULL))); if(_other->me->sig_upper != NULL && _other->me->sig_lower != NULL) { // ensure that we don't try and use the next node as something we will remove now if(_other->me->sig_upper->node->next == _other->me->sig_lower->node){ other = _other->me->sig_lower->node->next; } else { other = _other->me->sig_upper->node->next; } _delete_frequency_from_list(_other->me->sig_upper->node); _delete_frequency_from_list(_other->me->sig_lower->node); } } _delete_frequency_from_list(_other); } else other = other->next; } } curr = curr->next; } } void add_user_frequency(const char *name, double f){ assert(name != NULL); assert(strlen(name) < LINE_LEN); frequency_t *frequency = &inter.tmp_f_list[inter.num_tmp_frequencies]; frequency->f = f; frequency->index = inter.num_tmp_frequencies; frequency->type = USER_FIELD; frequency->carrier = NULL; frequency->order = 0; strcpy(frequency->name, name); frequency->isTuned = 0; inter.num_tmp_frequencies++; // increase number of frequencies _insert_frequency_to_list(NULL, frequency); } void add_carrier_frequency(const char *name, double f, frequency_t **freq){ assert(name != NULL); assert(strlen(name) < LINE_LEN); frequency_t *frequency = &inter.tmp_f_list[inter.num_tmp_frequencies]; frequency->f = f; frequency->index = inter.num_tmp_frequencies; frequency->type = CARRIER_FIELD; frequency->carrier = NULL; frequency->order = 0; strcpy(frequency->name, name); frequency->isTuned = is_put_tuned(&frequency->f, &frequency->tuning_component, &frequency->tuning_type); if(!frequency->isTuned) frequency->isTuned = is_xaxis_tuned(&frequency->f, &frequency->tuning_component, &frequency->tuning_type); if(!frequency->isTuned) { if(inter.minimizer != NULL && inter.minimizer->variable == &frequency->f) frequency->isTuned = true; } inter.num_tmp_frequencies++; // increase number of frequencies if(freq) *freq = frequency; _insert_frequency_to_list(NULL, frequency); } void add_sideband_frequency(const char *name, frequency_t **created_freq, frequency_t *car_freq, double *fmod, int type, int order){ assert(name != NULL); assert(strlen(name) < LINE_LEN); assert(type == MODULATOR_FIELD || type == SIGNAL_FIELD); assert(order != 0); // this is a sideband so must be something non zero assert(fmod != NULL); frequency_t *frequency = &inter.tmp_f_list[inter.num_tmp_frequencies]; frequency->f = car_freq->f + order * (*fmod); frequency->index = inter.num_tmp_frequencies; frequency->type = type; frequency->carrier = car_freq; frequency->order = order; if(type == SIGNAL_FIELD){ if(order == 1) car_freq->sig_upper = frequency; else if(order == -1) car_freq->sig_lower = frequency; else bug_error("unhandled signal sideband order"); } // if the carrier frequency is being tuned so will this sideband as the carrier // field will be changed if(car_freq->isTuned){ frequency->isTuned = car_freq->isTuned; frequency->tuning_component = NULL; frequency->tuning_type = CARRIER_FIELD; } // we arent't being tuned by the carrier our modulation frequency might be getting tuned if(!frequency->isTuned) frequency->isTuned = is_put_tuned(fmod, &frequency->tuning_component, &frequency->tuning_type); if(!frequency->isTuned) frequency->isTuned = is_put_tuned(&frequency->f, &frequency->tuning_component, &frequency->tuning_type); if(!frequency->isTuned) frequency->isTuned = is_xaxis_tuned(fmod, &frequency->tuning_component, &frequency->tuning_type); if(!frequency->isTuned) frequency->isTuned = is_xaxis_tuned(&frequency->f, &frequency->tuning_component, &frequency->tuning_type); if(!frequency->isTuned) frequency->isTuned = is_minimizer_tuned(&frequency->f, &frequency->tuning_component, &frequency->tuning_type); if(!frequency->isTuned) frequency->isTuned = is_minimizer_tuned(fmod, &frequency->tuning_component, &frequency->tuning_type); frequency->f_mod = fmod; strcpy(frequency->name, name); inter.num_tmp_frequencies++; // increase number of frequencies if(created_freq) *created_freq = frequency; frequency_node_t *cnode; if(!_find_frequency_in_list(frequency->carrier, &cnode)) bug_error("Could not find carrier in the frequency linked list"); _insert_frequency_to_list(cnode, frequency); } void __set_freq_cpling(char* name, double f_mod, int order, bool single_sb, bool numerical_f_couple, frequency_t *fin, frequency_t *fout, bool allocing, int ***f_couple_allocd, int ***f_couple_order, int ***does_f_couple, int i, int j, int k, double df){ // if either of the frequencies are being tuned we might // at some point couple between them if(allocing && (fin->isTuned || fout->isTuned)) f_couple_allocd[k][i][j] = 1; // check if this is a sideband and that its carrier is // the incoming frequency fin we are currently considering. // Then finally check that the modulation frequency of this // modulator is the same as the sidebands modulation frequency if(fout->carrier != NULL && fout->carrier == fin && f_mod == *fout->f_mod){ // if we do single sideband we only consider the positive sideband if(!single_sb || (single_sb && fout->order == 1)){ // fout is a sideband of fin if(allocing) f_couple_allocd[k][i][j] = 1; does_f_couple[k][i][j] = 1; f_couple_order[k][i][j] = fout->order; } } else if(numerical_f_couple) { // fout might not be a sb of fin but it could still couple int l; for(l=-order; l <= order; l++){ if(single_sb == true && l < 0) continue; if(df == l * f_mod){ // if both are static frequencies then see if we are harmonics // of the modulation frequency if(allocing) f_couple_allocd[k][i][j] = 1; does_f_couple[k][i][j] = 1; f_couple_order[k][i][j] = l; } if(!allocing && does_f_couple[k][i][j] && !f_couple_allocd[k][i][j]){ warn("enabled coupling between frequencies %i and %i at %s (modulation %.5gHz) but this coupling was not allocated for\n", i, j, name, f_mod); does_f_couple[k][i][j] = 0; f_couple_order[k][i][j] = 0; } } } } void __set_in_freq_cpling(int type, frequency_t *fin, frequency_t **f_list, bool allocing, int ***f_couple_allocd, int ***f_couple_order, int ***does_f_couple, int num, int num_freqs, int i){ assert(fin != NULL); assert(f_list != NULL); // Sets the coupling matrices to always couple from one frequency back into itself int j, k; for(k=0; k < num; k++){ if(allocing) f_couple_allocd[k][i][i] = 1; f_couple_order[k][i][i] = 0; does_f_couple[k][i][i] = 1; } for(j=0; jf - fin->f; // check all possible couplings at a modulator // ensure that carrier and modulator fields don't couple to signal fields as // they would need to be rescaled. However signal fields can couple between one // another if((fin->type != SIGNAL_FIELD && fout->type != SIGNAL_FIELD) || (fin->type == SIGNAL_FIELD && fout->type == SIGNAL_FIELD)){ if(type == MODULATOR) { for(k=0; k < inter.num_modulators; k++){ modulator_t *mod = &inter.modulator_list[k]; __set_freq_cpling(mod->name, mod->f_mod, mod->order, mod->single_sideband_mode, mod->numerical_f_couple, fin, fout, allocing, f_couple_allocd, f_couple_order, does_f_couple, i, j, k,df); } } else if(type == MIRROR) { for(k=0; k < inter.num_mirror_dithers; k++){ mirror_t *m = inter.mirror_dither_list[k]; __set_freq_cpling(m->name, m->dither_f, m->dither_order, false, true, fin, fout, allocing, f_couple_allocd, f_couple_order, does_f_couple, i, j, k, df); } } else if(type == BEAMSPLITTER) { for(k=0; k < inter.num_bs_dithers; k++){ beamsplitter_t *m = inter.bs_dither_list[k]; __set_freq_cpling(m->name, m->dither_f, m->dither_order, false, true, fin, fout, allocing, f_couple_allocd, f_couple_order, does_f_couple, i, j, k, df); } } else { bug_error("Component type not expected\n"); } } } } /** * loops through components that couple between frequencies and updates which * frequencies couple to which. This should be called if any of the frequencies are * changed. * * Set allocing to > 0 if this is being used to determine which frequencies will * couple for matrix and memory allocation */ void update_frequency_coupling(int allocing){ int g,i; ifo_matrix_vars_t *matrices[2] = {&M_ifo_car, &M_ifo_sig}; for (g=0; g<2; g++){ ifo_matrix_vars_t *matrix = matrices[g]; frequency_t **f_list = matrix->frequencies; int f_num = matrix->num_frequencies; for(i=0; i < f_num; i++){ frequency_t *fin = f_list[i]; assert(fin != NULL); // set all the direct frequency couplings for components __set_in_freq_cpling(MODULATOR, fin, f_list, allocing, matrix->mod_f_couple_allocd, matrix->mod_f_couple_order, matrix->mod_does_f_couple, inter.num_modulators, matrix->num_frequencies, i); __set_in_freq_cpling(MIRROR, fin, f_list, allocing, matrix->m_f_couple_allocd, matrix->m_f_couple_order, matrix->m_does_f_couple, inter.num_mirror_dithers, matrix->num_frequencies, i); __set_in_freq_cpling(BEAMSPLITTER, fin, f_list, allocing, matrix->bs_f_couple_allocd, matrix->bs_f_couple_order, matrix->bs_does_f_couple, inter.num_bs_dithers, matrix->num_frequencies, i); } } } void print_frequency_couplings(const char *name, int num_fs, int **does_f_couple, int **f_couple_order){ int i,j; message(" * Frequency coupling matrix for %s\n", name); message(" "); for(i=0; i < num_fs; i++) message("%-3i", i); message("\n"); for(i=0; i < num_fs; i++){ message("%-3i", i); for(j=0; j < num_fs; j++){ if(does_f_couple[i][j]){ message("%-3i", f_couple_order[i][j]); } else { message(" "); } } message("\n"); } } /** * Loops through each frequency, checks if its modulation or carrier's frequency * has been altered and if so updates itself */ void update_frequency_list(){ int i; if(carrier_frequencies_tuned){ for(i=0; itype){ case MODULATOR_FIELD: if(f->isTuned){ // If the sideband is tuned then the modulator frequency or // the carrier frequency is being updated. if(f->f_mod == NULL) bug_error("sideband did not have modulation frequency available for updating"); f->f = f->carrier->f + f->order * (*(f->f_mod)); } break; case SIGNAL_FIELD: bug_error("There shouldn't be any signal frequency in carrier vector"); case CARRIER_FIELD: case USER_FIELD: // if this is a carrier then nothing needs to be done as the // frequency pointer will have already been used to update its // value from puts or xaxis break; default: bug_error("unhandled signal field type"); } } } if(signal_frequencies_tuned){ for(i=0; itype){ case SIGNAL_FIELD: if(f->isTuned){ // If the sideband is tuned then the modulator frequency or // the carrier frequency is being updated. if(f->f_mod == NULL) bug_error("sideband did not have modulation frequency available for updating"); f->f = f->carrier->f + f->order * (*(f->f_mod)); } break; case MODULATOR_FIELD: case CARRIER_FIELD: bug_error("There shouldn't be any carrier or modulators in signal frequency vector"); case USER_FIELD: // if this is a carrier then nothing needs to be done as the // frequency pointer will have already been used to update its // value from puts or xaxis break; default: bug_error("unhandled signal field type"); } } } } void __fill_mod_freqs(frequency_t *curr_carrier, int single_sb, int order, char *name, char *fname, double *mod_freq){ if(single_sb == ON) assert(order == 1); int j, l, k; for(l=order; l > 0; l--){ for(k=-1; k <= 1; k+=2){ j = l * k; if(single_sb == ON && j <= 0) continue; if(j != 0){ sprintf(fname, "%s_%s_%s%i", name, curr_carrier->name, (j < 0) ? "m" : "p", (j<0)?-j:j); frequency_t *fmod = NULL; add_sideband_frequency(fname, &fmod, curr_carrier, mod_freq, MODULATOR_FIELD, j); if(inter.num_signals){ char _fname[LINE_LEN+20] = {0}; sprintf(_fname, "%s_msb", fmod->name); add_sideband_frequency(_fname, NULL, fmod, &inter.fsig, SIGNAL_FIELD, -1); sprintf(_fname, "%s_psb", fmod->name); add_sideband_frequency(_fname, NULL, fmod, &inter.fsig, SIGNAL_FIELD, 1); } } } } } /** * Loops through all fsig'd components, modulators and dither to add the frequencies * they generate */ void fill_frequency_linked_list(){ frequency_node_t *curr_carrier = inter.frequency_linked_list.head; frequency_node_t *next_carrier = NULL; // loop through every carrier field and add the sidebands while(curr_carrier != NULL){ int i; char fname[LINE_LEN + 10] = {0}; if(curr_carrier->me->type == USER_FIELD){ // add signal sidebands around manual frequency if(inter.num_signals){ sprintf(fname, "%s_msb", curr_carrier->me->name); add_sideband_frequency(fname, NULL, curr_carrier->me, &inter.fsig, SIGNAL_FIELD, -1); sprintf(fname, "%s_psb", curr_carrier->me->name); add_sideband_frequency(fname, NULL, curr_carrier->me, &inter.fsig, SIGNAL_FIELD, 1); } curr_carrier = curr_carrier->next; continue; } else if(curr_carrier->me->type != CARRIER_FIELD){ curr_carrier = curr_carrier->next; continue; } // save the next frequency now as we will add a bunch of sidebands next // after this carrier and there isn't much point iterating over them next_carrier = curr_carrier->next; // for every modulator add the sideband frequencies around this carrier for(i=0; ime, mod->single_sideband_mode, mod->order, mod->name, fname, &mod->f_mod); } for(i=0; ime, false, m->dither_order, m->name, fname, &m->dither_f); } for(i=0; ime, false, m->dither_order, m->name, fname, &m->dither_f); } // add in the signal sidebands if there are some fsig components if(inter.num_signals){ sprintf(fname, "%s_msb", curr_carrier->me->name); add_sideband_frequency(fname, NULL, curr_carrier->me, &inter.fsig, SIGNAL_FIELD, -1); sprintf(fname, "%s_psb", curr_carrier->me->name); add_sideband_frequency(fname, NULL, curr_carrier->me, &inter.fsig, SIGNAL_FIELD, 1); } curr_carrier = next_carrier; } } //! Get the signal type as a string /*! * The output will be one of SIG_AMP, SIG_PHS, SIG_FRQ, SIG_X or SIG_Y. * This can also be the input, but as variables. What is returned is a * stringified version of the input type's name. * * \param signal_type The signal type macro variable, or the type member of * a signal structure. * * \return The signal type as a string * * \see Test_get_signal_type_str() */ char *get_signal_type_str(int signal_type) { // check the input assert(signal_type == SIG_AMP || signal_type == SIG_PHS || signal_type == SIG_FRQ || signal_type == SIG_X || signal_type == SIG_Y); // generate the return string switch (signal_type) { case (SIG_AMP): return duplicate_string("SIG_AMP"); break; case (SIG_PHS): return duplicate_string("SIG_PHS"); break; case (SIG_FRQ): return duplicate_string("SIG_FRQ"); break; case (SIG_X): return duplicate_string("SIG_X"); break; case (SIG_Y): return duplicate_string("SIG_Y"); break; case (SIG_Z): return duplicate_string("SIG_Z"); case (SIG_SURF): return duplicate_string("SIG_SURF"); case (SIG_FZ): return duplicate_string("SIG_FZ"); case (SIG_FRX): return duplicate_string("SIG_FRX"); case (SIG_FRY): return duplicate_string("SIG_FRY"); case (SIG_FSURF): return duplicate_string("SIG_FSURF"); default: bug_error("Incorrect signal type"); } // return the empty string as the last-ditch option return duplicate_string(""); } /*! * * \return The signal units as a string * */ char *get_signal_unit_str(int component_type, int signal_type) { switch (signal_type) { case (SIG_AMP): switch (component_type){ case (MIRROR): case (BEAMSPLITTER): case (LIGHT_INPUT): return duplicate_string("sqrt(W)"); break; default: return duplicate_string("?"); break; } break; case (SIG_PHS): switch (component_type){ case (MIRROR): case (BEAMSPLITTER): case (LIGHT_INPUT): return duplicate_string("rad"); break; default: return duplicate_string("?"); break; } break; case (SIG_FRQ): return duplicate_string("Hz"); break; case (SIG_X): return duplicate_string("rad"); break; case (SIG_Y): return duplicate_string("rad"); break; case (SIG_Z): return duplicate_string("?"); case (SIG_SURF): return duplicate_string("?"); case (SIG_FZ): return duplicate_string("?"); case (SIG_FRX): return duplicate_string("?"); case (SIG_FRY): return duplicate_string("?"); case (SIG_FSURF): return duplicate_string("?"); default: bug_error("Incorrect signal type"); } // return the empty string as the last-ditch option return duplicate_string(""); } /* if (string_matches_exactly(signal_type_string, "phase")) { return SIG_PHS; } else if (string_matches_exactly(signal_type_string, "z")) { return SIG_Z; } else if (string_matches_exactly(signal_type_string, "phi")) { return SIG_PHS; } else if (string_matches_exactly(signal_type_string, "phs")) { return SIG_PHS; } else if (string_matches_exactly(signal_type_string, "amp")) { return SIG_AMP; } else if (string_matches_exactly(signal_type_string, "freq")) { return SIG_FRQ; } else if (string_matches_exactly(signal_type_string, "frq")) { return SIG_FRQ; } else if (string_matches_exactly(signal_type_string, "xbeta") || string_matches_exactly(signal_type_string, "rx") || string_matches_exactly(signal_type_string, "yaw")) { return SIG_X; } else if (string_matches_exactly(signal_type_string, "ybeta") || string_matches_exactly(signal_type_string, "ry") || string_matches_exactly(signal_type_string, "pitch")) { return SIG_Y; } else if (string_matches_exactly(signal_type_string, "p")) { return SIG_PHS; } else if (string_matches_exactly(signal_type_string, "a")) { return SIG_AMP; } else if (string_matches_exactly(signal_type_string, "f")) { return SIG_FRQ; } else if (string_matches_exactly(signal_type_string, "x")) { return SIG_X; } else if (string_matches_exactly(signal_type_string, "y")) { return SIG_Y; } else if (signal_type_string[0] == 's' && isdigit(signal_type_string[1])) { // check if we are trying to drive a surface mode return SIG_SURF; } 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; } else if ((signal_type_string[0] == 'F' && signal_type_string[1] == 's') && isdigit(signal_type_string[2])) { // check if we are trying to drive a surface mode return SIG_FSURF; } return (NOT_FOUND); } */ //! Check for too long line /*! * \param s input line to process * * \see Test_prepare_line() */ void check_for_long_line(FILE *fp, char *s) { char *sc; // don't do test at the end of the file if (!feof(fp)){ int found_end_of_line = 0; if ((sc = strchr(s, '\n')) != NULL) { found_end_of_line = 1; } if ((sc = strchr(s, '\r')) != NULL) { found_end_of_line = 1; } if (NOT found_end_of_line) { warn("line '%.20s ...': did not find end of line, maybe line to long (>512 chars)?\n", s); } } // check that nothing in the line is longer than MAX_TOKEN_LEN-1 chararcters char tmp_line[LINE_LEN]; strcpy(tmp_line, s); char *tmp_pointer; char *ap; int max_length = 0; int new_length; tmp_pointer = tmp_line; ap = strtok(tmp_pointer, " \t"); while (ap != NULL) { if (*ap != '\0') { new_length = strlen(ap); if (new_length > max_length) { max_length = new_length; } } ap = strtok(NULL, " \t"); } if (max_length >= MAX_TOKEN_LEN) { gerror("line '%.20s ...': a part of this line is larger than %d characters.\n", s, MAX_TOKEN_LEN - 1); } } //! Prepare input line for processing /*! * \param s input line to process * \param mode if mode == 1, treat " as a comment character * * \see Test_prepare_line() */ int prepare_line(char *s, int mode) { // mode=1 removes '"' as comments char *sc; // sanity checks on inputs // mode should be either 0 or 1 assert(mode == 0 || mode == 1); if ((sc = strchr(s, '\n')) != NULL) { *sc = '\0'; /* remove linefeed */ } if ((sc = strchr(s, '\r')) != NULL) { *sc = '\0'; /* remove carriage return */ } if (s[0] == '#' || s[0] == '%' || s[0] == '"') { return (1); } if (s[0] == '\n' || s[0] == '\r' || s[0] == '\0') { return (1); } sc = strchr(s, '#'); if (sc != NULL) { *sc = '\0'; /* remove comment from end of line */ } sc = strchr(s, '%'); if (sc != NULL) { *sc = '\0'; /* remove comment from end of line */ } if (mode) { sc = strchr(s, '"'); if (sc != NULL) { *sc = '\0'; // remove comment from end of line } } strcat(s, " "); /* add space at end for sscanf(...%n) */ return (0); } //! Check the function name /*! * \param function_name the function name to check * \param s the command string (i.e. the entire line) from the kat file * * \todo should this be in kat_check.c ??? * * \see Test_check_funcname() */ void check_funcname(char *function_name, const char *input_string) { if ((strcasecmp(function_name, "x") == 0) || (strcasecmp(function_name, "x1") == 0) || (strcasecmp(function_name, "x2") == 0) || (strcasecmp(function_name, "x3") == 0)) { gerror("line '%s':\nx/x1/x2/x3 not allowed as set names.\n", input_string); } bool cmd_name_is_found = false; int func_cmd_index; for (func_cmd_index = 0; func_cmd_index < inter.num_func_cmds; func_cmd_index++) { if (strcasecmp(function_name, inter.function_list[func_cmd_index].name + 1) == 0) { cmd_name_is_found = true; } } if (cmd_name_is_found) { gerror("name '%s' already used for a function.\n", input_string); } int set_cmd_index; cmd_name_is_found = false; for (set_cmd_index = 0; set_cmd_index < inter.num_set_cmds; set_cmd_index++) { if (strcasecmp(function_name, inter.set_list[set_cmd_index].name + 1) == 0) { cmd_name_is_found = true; } } if (cmd_name_is_found) { gerror("name '%s' already used for a \'set\'.\n", input_string); } cmd_name_is_found = false; int lock_index; for (lock_index = 0; lock_index < inter.num_locks; lock_index++) { if (strcasecmp(function_name, inter.lock_list[lock_index].name + 1) == 0) { cmd_name_is_found = true; } } if (cmd_name_is_found) { gerror("name '%s' already used for a lock.\n", input_string); } } //! Find the equals sign /*! * \param function_name gets the part of the string before the equals sign * \param function_string the string to search * * \see Test_find_equal() */ int find_equal(char *function_name, char *function_string) { int error_code = 0; char *equals_sign_pointer; char temporary_function_name[LINE_LEN] = {0}; strcpy(temporary_function_name, function_string); equals_sign_pointer = strchr(temporary_function_name, '='); if (equals_sign_pointer != NULL) { // copy second part of string into function_name and remove it from // original part of string strcpy(function_string, equals_sign_pointer + 1); *equals_sign_pointer = '\0'; if (strlen(temporary_function_name) > LINE_LEN) { error_code = 2; } else { strcpy(function_name, temporary_function_name); } } else { error_code = 1; } return (error_code); } //! Sort variables so that the longer ones come first /*! * (helps to avoid problems with similar names like "var1" "var12") * * \see Test_sort_variable_names() */ void sort_variable_names(void) { constant_t vtmp; int i, j; for (i = 0; i < inter.num_constants; i++) { for (j = inter.num_constants; j > i; j--) { if (strlen(inter.constant_list[i].name) < strlen(inter.constant_list[j].name)) { vtmp = inter.constant_list[i]; inter.constant_list[i] = inter.constant_list[j]; inter.constant_list[j] = vtmp; } } } } //! Sort set names so that the longer ones come first /*! * (helps to avoid problems with similar names like "var1" "var12") * * \todo untested * \todo sort_set -> sort_set_cmd_names * * \see Test_sort_set() */ void sort_set(void) { set_command_t xtmp; int i, j; for (i = 0; i < inter.num_set_cmds; i++) { for (j = inter.num_set_cmds; j > i; j--) { if (strlen(inter.set_list[i].name) < strlen(inter.set_list[j].name)) { xtmp = inter.set_list[i]; inter.set_list[i] = inter.set_list[j]; inter.set_list[j] = xtmp; } } } } //! Replace variable names by values /*! * \param input_string input line read in by Finesse * * \see Test_insert_constants_values() */ void insert_constants_values(char *input_string) { int constant_index; for (constant_index = 0; constant_index < inter.num_constants; constant_index++) { int num_strings_replaced = replace_string( input_string, inter.constant_list[constant_index].name, inter.constant_list[constant_index].value, inter.constant_list[constant_index].name_len, inter.constant_list[constant_index].value_len); if (num_strings_replaced == -1) { gerror("line '%s':\nreplacing constants exceeds line-length limit\n", input_string); } } } //! Search for s1 in s and replace by s2 /*! * Check if resulting string is larger then LINE_LEN * return -1 on error or number of replaced strings * * \param input_string input line * \param search_string the string to search for * \param replacement_string the replacement string * \param search_string_len the length of s1 * \param replacement_string_len the length of s2 * * \todo s0, sc, sct, st -> better names * * \see Test_replace_string() */ int replace_string(char *input_string, char *search_string, char *replacement_string, int search_string_len, int replacement_string_len) { assert(input_string != NULL); assert(search_string != NULL); assert(replacement_string != NULL); assert(strlen(replacement_string) == replacement_string_len); //assert(strlen(search_string) == search_string_len); size_t in_len = strlen(input_string); char *input_string_copy = (char*) calloc(sizeof(char), in_len+replacement_string_len); char *st = (char*) calloc(sizeof(char), in_len+replacement_string_len); char *str_tmp = (char*) calloc(sizeof(char),in_len+replacement_string_len); if(input_string_copy == NULL || st == NULL || str_tmp == NULL) { gerror("Memory allocation error in replace_string\n"); } char *sc = NULL, *sct = NULL; int num_strings_replaced = 0; int offset_len = 0; strcpy(input_string_copy, input_string); sc = input_string_copy; while (sc != NULL) { sc = strstr(input_string_copy, search_string); if (sc != NULL) { num_strings_replaced++; offset_len += replacement_string_len - search_string_len; strcpy(str_tmp, sc + search_string_len); strcpy(input_string_copy, str_tmp); } } if (num_strings_replaced > 0) { strcpy(input_string_copy, input_string); strcpy(input_string, ""); sc = input_string_copy; while (sc != NULL) { sct = input_string_copy; sc = strstr(input_string_copy, search_string); if (sc != NULL) { // s to st 130905 strncpy(st, sct, strlen(sct) - strlen(sc)); st[strlen(sct) - strlen(sc)] = '\0'; strcat(input_string, st); strcat(input_string, replacement_string); strcpy(str_tmp, sc + search_string_len); strcpy(input_string_copy, str_tmp); } else { if (input_string != sct) { strcat(input_string, sct); } } } } free(input_string_copy); free(st); free(str_tmp); return (num_strings_replaced); } //! For given map name return the index /*! * \param mirror_index index of the mirror in which the map is to be found * \param mapname string containing the map name * */ int get_mirror_index_from_name(int mirror_index, char *mapname) { int i; for (i = 0; i < inter.mirror_list[mirror_index].num_maps; i++) { if (strcmp(mapname, inter.mirror_list[mirror_index].map[i]->name) == 0) { return i; } } return NOT_FOUND; } //! For given n, m return field index /*! * \param n mode of TEM field * \param m mode of TEM field * * \see Test_get_field_index_from_tem() * \see kat_inline for the reverse function * TODO: maybe check for inter.tem_is_set and either throw error or return field=0 for n=m=0 and field=1 otherwise? */ int get_field_index_from_tem(int n, int m) { int field; // sanity checks on inputs // n and m both must be greater than or equal to zero assert(m >= 0); assert(n >= 0); field = n * (inter.tem + 2) - (n * (n + 1)) / 2 + m; return (field); } //! For given n, m return field index for LG beam /*! * \param p mode of TEM field * \param l mode of TEM field */ int get_field_index_from_tem_LG(int p, int l) { assert(l <= inter.tem); assert(l >= -inter.tem); assert(p <= inter.tem); assert(p >= 0); int O = 2*p + abs(l); return (O*(O+2)+l)/2; } //! For given n, m and a maximum order maxtem, return field index /*! * \param n mode of TEM field * \param m mode of TEM field * */ int get_field_index_from_local_tem(int n, int m, int maxtem) { int field; // sanity checks on inputs // n, m and maxtem must be greather than or equal to zero assert(m >= 0); assert(n >= 0); assert(maxtem >= 0); field = n * (maxtem + 2) - (n * (n + 1)) / 2 + m; return (field); } //! For given field number return n, m as in TEM_nm /*! * * * \see Test_set_all_tems() */ void set_all_tems(void) { int k; assert(mem.all_tem_HG != NULL); assert(mem.all_tem_LG != NULL); int n, m; for (n = 0; n <= inter.tem; n++) { for (m = 0; m <= (inter.tem - n); m++) { if ((n + m) <= inter.tem) { k = get_field_index_from_tem(n, m); mem.all_tem_HG[2 * k] = n; mem.all_tem_HG[2 * k + 1] = m; } } } int p,l; for (p=0; p<=inter.tem; p++){ for (l=-inter.tem; l<=inter.tem; l++){ if (2*p + abs(l)<= inter.tem){ k = get_field_index_from_tem_LG(p, l); mem.all_tem_LG[2 * k] = p; mem.all_tem_LG[2 * k + 1] = l; } } } } /*** * Gets the complex beam parameter at a node going into a component. * Doesn't do any checks to see if the node is connected to the component. * * node: ptr to node struct * component_idx: unique component index * qx, qy: ptr to complex_t to store beam params in, if NULL it won't set * * Returns 0 if ground/dump node */ int get_node_q_to(node_t *node, int component_idx, complex_t *qx, complex_t *qy){ if (NOT node->gnd_node) { // If q was traced away from a component, a node it will have the component index // set to it, so here we reverse it if (node->component_index == component_idx) { if(qx) *qx = cminus(cconj(node->qx)); if(qy) *qy = cminus(cconj(node->qy)); } else { if(qx) *qx = node->qx; if(qy) *qy = node->qy; } } else { return 0; } } /*** * Gets the complex beam parameter at a node going from a component. * Doesn't do any checks to see if the node is connected to the component. * * node: ptr to node struct * component_idx: unique component index * qx, qy: ptr to complex_t to store beam params in, if NULL it won't set * * Returns 0 if ground/dump node */ int get_node_q_from(node_t *node, int component_idx, complex_t *qx, complex_t *qy){ if (NOT node->gnd_node) { // If q was traced "from" a node it will have the component index // set to it if (node->component_index == component_idx) { if(qx) *qx = node->qx; if(qy) *qy = node->qy; } else { if(qx) *qx = cminus(cconj(node->qx)); if(qy) *qy = cminus(cconj(node->qy)); } } else { return 0; } } //! Returns factorial for n<=100 /*! * using the fact that \f$\log(n!) = \log(n) + \log((n+1)!)\f$ * the \f$\sum_n \log(n)\f$ are stored in a static list * (natLogsOfN) while they are build. * * \param n order of factorial to find * * \see Test_fac() */ double fac(int n) { // sanity check on input assert(n >= 0); if (n > 100) { bug_error("factorial index too large"); } if (n <= natLogMax) { return (exp(natLogsOfN[n])); } int i; for (i = natLogMax + 1; i <= n; i++) { natLogsOfN[i] = natLogsOfN[i - 1] + log(i); } natLogMax = n; return (exp(natLogsOfN[n])); } //! Prints the (elapsed) time given in seconds in a nicer format /*! * \param elapsed_seconds the time to print * * \see Test_print_time() */ void print_time() { double elapsed_seconds = 0.0; timespec_t end; getTime(&end); //#if OS == __WIN_32__ // elapsed_seconds = end.tv_sec - clock_start.tv_sec + (end.tv_usec - clock_start.tv_usec) / 1e6; //#else elapsed_seconds = end.tv_sec - clock_start.tv_sec + (end.tv_nsec - clock_start.tv_nsec) / 1e9; //#endif char message_str[ERR_MSG_LEN]; message("computation time: "); sprint_time(message_str, elapsed_seconds); message("%s\n", message_str); } //! returns the (elapsed) time given in seconds in a string of `nicer' format /*! * \param elapsed_seconds the time to print * * \see Test_print_time() */ void sprint_time(char *message_str, double elapsed_seconds) { double integer_part, fractional_part; int hours = 0; int minutes = 0; int seconds = 0; if (options.perl1) { sprintf(message_str, "%.6gs", elapsed_seconds); } else if (elapsed_seconds > 3600) { fractional_part = modf(elapsed_seconds / 3600.0, &integer_part); hours = (int) integer_part; fractional_part = modf((elapsed_seconds / 60. - hours * 60), &integer_part); minutes = (int) integer_part; seconds = floor(elapsed_seconds - hours * 3600 - minutes * 60 + 0.5); sprintf(message_str, "%dh%dm%ds", hours, minutes, seconds); } else if (elapsed_seconds > 60) { fractional_part = modf(elapsed_seconds / 60.0, &integer_part); minutes = (int) integer_part; seconds = floor(elapsed_seconds - minutes * 60 + 0.5); sprintf(message_str, "%dm%ds", minutes, seconds); } else { if (elapsed_seconds < 1) { sprintf(message_str, "%.6gs", elapsed_seconds); } else { sprintf(message_str, "%.3gs", elapsed_seconds); } } (void) fractional_part; // silence compiler warning. } //! Returns input in decibels /*! * \param x input value to convert to decibels * * \see Test_xdb() */ double xdb(double x) { return 20 * log10(x); } //! Hermite Polynomial; Siegmann S.686 /*! * This function calculates the Hermite polynomial using the recursion * relation: * \f[ * H_{n+1}(x) = 2x H_n(x) - 2n H_{n-1}(x) * \f] * where the first four terms are given by: * \f{eqnarray*} * H_0 &=& 1\\ * H_1(x) &=& 2x\\ * H_2(x) &=& 4x^2 - 2\\ * H_3(x) &=& 8x^3 - 12x * \f} * * \param n the order of the Hermite polynomial * \param x position on x axis * * \todo not completely tested * * \see Test_hermite() */ double hermite(int n, double x) { // sanity check on input // n must be greater than or equal to zero assert(n >= 0); switch (n) { case 0: return (1.0); break; case 1: return (2.0 * x); break; case 2: return (4.0 * x * x - 2.0); break; case 3: return (8.0 * x * x * x - 12.0 * x); break; case 4: return (16.0 * gsl_pow_4(x) - 48.0 * x * x + 12.0); break; case 5: return (32.0 * gsl_pow_5(x) - 160.0 * x * x * x + 120.0 * x); break; case 6: return (64.0 * gsl_pow_6(x) - 480.0 * gsl_pow_4(x) + 720.0 * x * x - 120.0); break; case 7: return (128.0 * gsl_pow_7(x) - 1344.0 * gsl_pow_5(x) + 3360.0 * x * x * x - 1680.0 * x); break; case 8: return (256.0 * gsl_pow_8(x) - 3584.0 * gsl_pow_6(x) + 13440.0 * gsl_pow_4(x) - 13440.0 * x * x + 1680.0); break; case 9: return (512.0 * gsl_pow_9(x) - 9216.0 * gsl_pow_7(x) + 48384.0 * gsl_pow_5(x) - 80640.0 * x * x * x + 30240.0 * x); break; case 10: return (1024.0 * gsl_pow_2(gsl_pow_5(x)) - 23040.0 * gsl_pow_8(x) + 161280.0 * gsl_pow_6(x) - 403200.0 * gsl_pow_4(x) + 302400.0 * x * x - 30240.0); break; case 11: return -665280 * x + 2217600 * gsl_pow_3(x) - 1774080 * gsl_pow_5(x) + 506880 * gsl_pow_7(x) - 56320 * gsl_pow_9(x) + 2048 * gsl_pow_2(gsl_pow_5(x))*x; case 12: return 4096 * gsl_pow_2(gsl_pow_6(x)) - 135168 * gsl_pow_2(gsl_pow_5(x)) + 1520640 * gsl_pow_8(x) - 7096320 * gsl_pow_6(x) + 13305600 * gsl_pow_4(x) - 7983360 * x*x + 665280; case 13: return 8192 * gsl_pow_2(gsl_pow_6(x))*x - 319488 * gsl_pow_2(gsl_pow_5(x))*x + 4392960 * gsl_pow_9(x) - 26357760 * gsl_pow_7(x) + 69189120 * gsl_pow_5(x) - 69189120 * gsl_pow_3(x) + 17297280 * x; case 14: return 16384 * gsl_pow_2(gsl_pow_7(x))-745472 * gsl_pow_2(gsl_pow_6(x))+12300288 * gsl_pow_2(gsl_pow_5(x))-92252160 * gsl_pow_8(x)+322882560 * gsl_pow_6(x)-484323840 * gsl_pow_4(x)+242161920 * x*x-17297280; case 15: return 32768 * x*gsl_pow_2(gsl_pow_7(x))-1720320 * x*gsl_pow_2(gsl_pow_6(x))+33546240 * x*gsl_pow_2(gsl_pow_5(x))-307507200 * gsl_pow_9(x)+1383782400 * gsl_pow_7(x)-2905943040 * gsl_pow_5(x)+2421619200 * gsl_pow_3(x)-518918400 * x; case 16: return 65536 * gsl_pow_2(gsl_pow_8(x))-3932160 * gsl_pow_2(gsl_pow_7(x))+89456640 * gsl_pow_2(gsl_pow_6(x))-984023040 * gsl_pow_2(gsl_pow_5(x))+5535129600 * gsl_pow_8(x)-15498362880 * gsl_pow_6(x)+19372953600 * gsl_pow_4(x)-8302694400 * x*x+518918400; case 17: return 131072 * x*gsl_pow_2(gsl_pow_8(x))-8912896 * x*gsl_pow_2(gsl_pow_7(x))+233963520 * x*gsl_pow_2(gsl_pow_6(x))-3041525760 * x*gsl_pow_2(gsl_pow_5(x))+20910489600 * gsl_pow_9(x)-75277762560 * gsl_pow_7(x)+131736084480 * gsl_pow_5(x)-94097203200 * gsl_pow_3(x)+17643225600 * x; case 18: return 262144 * gsl_pow_2(gsl_pow_9(x))-20054016 * gsl_pow_2(gsl_pow_8(x))+601620480 * gsl_pow_2(gsl_pow_7(x))-9124577280 * gsl_pow_2(gsl_pow_6(x))+75277762560 * gsl_pow_2(gsl_pow_5(x))-338749931520 * gsl_pow_8(x)+790416506880 * gsl_pow_6(x)-846874828800 * gsl_pow_4(x)+317578060800 * x*x-17643225600; case 19: return 524288 * x*gsl_pow_2(gsl_pow_9(x))-44826624 * x*gsl_pow_2(gsl_pow_8(x))+1524105216 * x*gsl_pow_2(gsl_pow_7(x))-26671841280 * x*gsl_pow_2(gsl_pow_6(x))+260050452480 * x*gsl_pow_2(gsl_pow_5(x))-1430277488640 * gsl_pow_9(x)+4290832465920 * gsl_pow_7(x)-6436248698880 * gsl_pow_5(x)+4022655436800 * gsl_pow_3(x)-670442572800 * x; case 20: return 1048576 * gsl_pow_4(gsl_pow_5(x))-99614720 * gsl_pow_2(gsl_pow_9(x))+3810263040 * gsl_pow_2(gsl_pow_8(x))-76205260800 * gsl_pow_2(gsl_pow_7(x))+866834841600 * gsl_pow_2(gsl_pow_6(x))-5721109954560 * gsl_pow_2(gsl_pow_5(x))+21454162329600 * gsl_pow_8(x)-42908324659200 * gsl_pow_6(x)+40226554368000 * gsl_pow_4(x)-13408851456000 * x*x+670442572800; default: return (2 * x * hermite(n - 1, x) - 2 * (n - 1) * hermite(n - 2, x)); } bug_error("hermite 1"); return 0; } //! Bessel function of the first kind /*! * \param n order of the Bessel function * \param x position on x-axis * * \see Test_bessn() */ double bessn(int n, double x) { int j, m; double x1, bj, bjm, bjp, sum, tox, a; // sanity check on inputs // n must be greater than or equal to zero assert(n >= 0); if (n < 0) { bug_error("index < 0"); } else if (n == 0) { a = dbesj0(x); return a; } else if (n == 1) { a = dbesj1(x); return a; } else { x1 = fabs(x); tox = 2.0 / x1; if (x1 == 0) { return 0.0; } else if (x1 > (float) n) { bjm = dbesj0(x1); bj = dbesj1(x1); int i; for (i = 1; i < n; i++) { bjp = i * tox * bj - bjm; bjm = bj; bj = bjp; } a = bj; } else { m = 2 * ((n + (int) sqrt(40 * n)) / 2); // 40 defines accuracy of result j = 0; bjp = a = sum = 0.0; bj = 1.0; int i; for (i = m; i > 0; i--) { bjm = i * tox * bj - bjp; bjp = bj; bj = bjm; if (fabs(bj) > 1.0e10) { bj *= 1.0e-10; bjp *= 1.0e-10; a *= 1.0e-10; sum *= 1.0e-10; } if (j) { sum += bj; } j = !j; if (i == n) { a = bjp; } } sum = 2.0 * sum - bj; a /= sum; } return (x < 0.0 && (n & 1) ? -a : a); } return 0.0; /* silence compiler warning */ } //! Check the input name of a component /*! * \param s input line from .kat file to check * * \see Test_check_name() */ void check_name(const char *line, const char *input_string) { int found = 0; double dtmp; // check if name is a number: this indicates that the name has been // forgotten and a first parameter is treated as a name // e.g. a space: s 10 1 n1 n2 if (!(atod((char *) input_string, &dtmp))) { gerror("Line %s:\nnumeric name '%s' not allowed\n", line, input_string); } int i; for (i = 0; i < num_names; i++) { if (strcasecmp(input_string, mem.names[i]) == 0) { found = 1; } } if (found) { warn("Line %s:\nWarning: name '%s' already used. This may cause confusion.\n", line, input_string); } else { assert(num_names < mem.num_node_names); mem.names[num_names++] = duplicate_string(input_string); } } //! Find a component name (and check that it's not used more than once) /*! * \param component_name_string the current line of the .kat file * * \todo function not used - remove??? * * \todo untested * * \see Test_find_name() */ int find_name(const char *component_name_string) { int num = -1; int found = 0; int name_index; for (name_index = 0; name_index < num_names; name_index++) { if (strcasecmp(component_name_string, mem.names[name_index]) == 0) { found++; num = name_index; } } if (found > 1) { gerror("component name '%s' used twice.\n", component_name_string); } return (num); } //! Duplicate string /*! * \param input_string input string to duplicate * * \return duplicated_string the duplicated string * * \see Test_duplicate_string() */ static int rest_string_length = 0; //!< length of the rest of the string static char *next_string = NULL; //!< the next string to consider static char *start_string = NULL; //!< the start string??? Not used, remove? char *duplicate_string(const char *input_string) { int string_length = strlen(input_string) + 1; char *duplicated_string; if (string_length >= ALLOC_LEN) { duplicated_string = (char *) malloc(string_length); if (duplicated_string == NULL) { gerror("malloc1\n"); } memcpy(duplicated_string, input_string, string_length); return duplicated_string; } if (string_length > rest_string_length) { /* if (start_string != NULL) dupdump(); */ start_string = (char *) malloc(ALLOC_LEN); next_string = (char *) malloc(ALLOC_LEN); if (next_string == NULL) { gerror("malloc2\n"); } rest_string_length = ALLOC_LEN; } duplicated_string = next_string; next_string += string_length; rest_string_length -= string_length; memcpy(duplicated_string, input_string, string_length); return duplicated_string; } //! Duplicate string and allocate n extra bytes /*! * \param input_string input string to duplicate * \param num_extra_bytes the number of extra bytes to allocate * * \return return_string * * \see Test_duplicate_to_longer_string() */ char *duplicate_to_longer_string(const char *input_string, int num_extra_bytes) { int string_length = strlen(input_string) + 1; int new_string_length = string_length + num_extra_bytes; char *return_string; // sanity check on input assert(num_extra_bytes >= 0); if (new_string_length >= ALLOC_LEN) { return_string = malloc(new_string_length); if (return_string == NULL) { gerror("malloc1a\n"); } memcpy(return_string, input_string, string_length); return return_string; } if (new_string_length > rest_string_length) { if (start_string != NULL) { start_string = next_string = malloc(ALLOC_LEN); } if (next_string == NULL) { gerror("malloc2b\n"); } rest_string_length = ALLOC_LEN; } return_string = next_string; next_string += new_string_length; rest_string_length -= new_string_length; memcpy(return_string, input_string, string_length); return return_string; } //! Given a node name return the node index /*! * \param node_name input string of node name * * \return the node index or NODE_UNUSED if the node is not used. * * \see Test_get_node_index_for() */ int get_node_index_for(const char *node_name) { /* existing node */ int node_index; for (node_index = 0; node_index < inter.num_nodes; node_index++) { if (strcmp(node_name, inter.node_list[node_index].name) == 0) { debug_msg("Known node found: %s\n", node_name); return node_index; } } /* otherwise it's not used */ return NODE_UNUSED; } const char *get_comp_name2(int type, int list_index){ switch(type){ case MIRROR: return inter.mirror_list[list_index].name; case BEAMSPLITTER: return inter.bs_list[list_index].name; case LENS: return inter.lens_list[list_index].name; case SAGNAC: return inter.sagnac_list[list_index].name; case MODULATOR: return inter.modulator_list[list_index].name; case DIODE: return inter.diode_list[list_index].name; case SPACE: return inter.space_list[list_index].name; default: bug_error("type not handled"); } return NULL; } const char *get_comp_name(int type, void *comp){ switch(type){ case MIRROR: return ((mirror_t*)comp)->name; case BEAMSPLITTER: return ((beamsplitter_t*)comp)->name; case LENS: return ((lens_t*)comp)->name; case SAGNAC: return ((sagnac_t*)comp)->name; case MODULATOR: return ((modulator_t*)comp)->name; case DIODE: return ((diode_t*)comp)->name; case SPACE: return ((space_t*)comp)->name; default: bug_error("type not handled"); } return NULL; } /* Stores the components that are attached to the node, only two can * be attached to each node. */ void connect_component_to_node(int node_idx, void *comp, int type){ assert(comp != NULL); assert(node_idx >= 0 && node_idx < mem.num_nodes); node_connections_t *nc = &inter.node_conn_list[node_idx]; assert(nc != NULL); if(nc->comp_1 == NULL) { nc->comp_1 = comp; nc->type_1 = type; } else if(nc->comp_1 == comp){ gerror("%s is already connected to node %s\n", get_comp_name(nc->type_1, nc->comp_1), inter.node_list[node_idx].name); } else if(nc->comp_2 == NULL && !(inter.node_list[node_idx].io & 2)) { // need check to ensure that there aren't more than 1 component at an input node nc->comp_2 = comp; nc->type_2 = type; } else gerror("Tried to attach more than 2 components to node %s\n", inter.node_list[node_idx].name); if(nc->comp_1 != NULL && nc->comp_2 != NULL){ if(nc->type_1 != SPACE && nc->type_2 != SPACE){ const char *name1 = get_comp_name(nc->type_1, nc->comp_1); const char *name2 = get_comp_name(nc->type_2, nc->comp_2); gerror("Components %s and %s must be separated by a space from version 1.1 of Finesse\n", name1, name2); } } } //! For a given node name update and return node index /*! * \param node_name input string of node name * * \return the updated (or new) node index * * \see Test_update_node_index_for() */ int update_node_index_for(const char *node_name, int direction) { assert(direction == IN_OUT || direction == OUT_IN || direction == EITHER_DIR); assert(node_name != NULL); /* existing node */ int node_index; for (node_index = 0; node_index < inter.num_nodes; node_index++) { if (strcmp(node_name, inter.node_list[node_index].name) == 0) { debug_msg("Known node found: %s\n", node_name); // if the direction hasn't been set preciously, i.e. a laser input node // and this component is trying to set the direction, then let it do so if(inter.node_list[node_index].direction == EITHER_DIR && direction != EITHER_DIR) inter.node_list[node_index].direction = direction; // if the component is willing to accept either direction, i.e and output node // then don't worry, otherwise throw error //if(direction != EITHER_DIR && inter.node_list[node_index].direction != direction) // bug_error("Tried to update node with different direction to what has already been set"); return node_index; } } /* if too many nodes, barf */ if (inter.num_nodes >= mem.num_nodes) { bug_error("update node index: too many nodes\n"); } /* else new node */ node_t *node; node = &inter.node_list[inter.num_nodes]; node->list_index = node_index; node->direction = direction; // a new dump node if (strcasecmp(node_name, "dump") == 0) { // make up a new name for the dump mode. // format: dump_xxxxx, where x is an integer char new_dump_name[20] = {0}; sprintf(new_dump_name, "dump_%05d", inter.num_dump_nodes); debug_msg("New dump node created: %s\n", new_dump_name); // increment the number of dumps in the interferometer inter.num_dump_nodes++; // check its name to make sure things don't clash check_name("Automatic dump note", new_dump_name); // assign the name to the node strcpy(node->name, new_dump_name); // initialise connection and io-node information node->connect = 0; node->io = 0; // initialise the index of refraction to the default for a new node! //node->n = &(init.n0); // initialse qset for new node, i.e. q value is not set yet! node->q_is_set = false; // initialse qx and qy to zero node->qx = complex_0; node->qy = complex_0; // it's a dump node so set the gnd_node attribute to true node->gnd_node = true; // calculate the number of reduced nodes inter.num_reduced_nodes = inter.num_nodes + 1 - inter.num_dump_nodes; // and return the node index return (inter.num_nodes++); }// a new normal node else { debug_msg("New normal node created: %s\n", node_name); check_name("Node name check", node_name); // give the node a name strcpy(node->name, node_name); // initialise connection and io-node information node->connect = 0; node->io = 0; // initialise the index of refraction to the default for a new node! node->n = &(init.n0); // initialse qset for new node, i.e. q value is not set yet! node->q_is_set = false; // initialse qx and qy to zero node->qx = complex_0; node->qy = complex_0; // it's not a gnd_node, so set that attribute to false node->gnd_node = false; // calculate the number of reduced nodes inter.num_reduced_nodes = inter.num_nodes + 1 - inter.num_dump_nodes; // and return the node index return (inter.num_nodes++); } } //! Calculate the reduced node index for the given node index /*! * To have a reduced (classical) matrix as a result of using dump nodes, one * needs to calculate what the correct node index would be were the dump * nodes not treated as "normal" nodes, this then gives the correct * calculation of the interferometer matrix indices. The reduced node * index is the same value as what occurred in older versions of Finesse * where one set dump nodes to have a node index of -1. * * This routine searches through the list of nodes and counts the number of * dump nodes before the given node index and returns the difference between * the node index and the number of dump nodes found. * * \param search_node_index The node index to search to. * * \return The reduced node index * * \see Test_get_reduced_index_for() */ int get_reduced_index_for(int search_node_index) { int num_dump_nodes_found = 0; // check the input assert(search_node_index >= 0); // if there aren't any dump nodes, return directly with the node_index if (inter.num_dump_nodes == 0) { return search_node_index; } else { // count up to the current node index in the node list and look for // the number of dumps before the current node index int node_index; for (node_index = 0; node_index <= search_node_index; node_index++) { if (inter.node_list[node_index].gnd_node) { num_dump_nodes_found++; } } return (search_node_index - num_dump_nodes_found); } } //! Given a component number with respect to the type, return overall number /*! * \param type number specifying the type of component * \param component_index the component index * * \see Test_get_overall_component_index() */ int get_overall_component_index(int type, int component_index) { // sanity checks on inputs assert(component_index >= 0); assert(component_index < inter.num_components); int overall_component_index = 0; switch (type) { case MIRROR: overall_component_index = 0; break; case BEAMSPLITTER: overall_component_index = inter.num_mirrors; break; case GRATING: overall_component_index = inter.num_mirrors + inter.num_beamsplitters; break; case SPACE: overall_component_index = inter.num_mirrors + inter.num_beamsplitters + inter.num_gratings; break; case LENS: overall_component_index = inter.num_mirrors + inter.num_beamsplitters + inter.num_gratings + inter.num_spaces; break; case DIODE: overall_component_index = inter.num_mirrors + inter.num_beamsplitters + inter.num_gratings + inter.num_spaces + inter.num_lenses; break; case LIGHT_INPUT: overall_component_index = inter.num_mirrors + inter.num_beamsplitters + inter.num_spaces + inter.num_lenses + inter.num_gratings + inter.num_diodes; break; case MODULATOR: overall_component_index = inter.num_mirrors + inter.num_beamsplitters + inter.num_spaces + inter.num_lenses + inter.num_gratings + inter.num_diodes + inter.num_light_inputs; break; case DBS: overall_component_index = inter.num_mirrors + inter.num_beamsplitters + inter.num_spaces + inter.num_lenses + inter.num_gratings + inter.num_diodes + inter.num_light_inputs + inter.num_modulators; break; case FSIG: overall_component_index = inter.num_mirrors + inter.num_beamsplitters + inter.num_spaces + inter.num_lenses + inter.num_gratings + inter.num_diodes + inter.num_light_inputs + inter.num_modulators + inter.num_dbss; break; case VARIABLE: overall_component_index = inter.num_mirrors + inter.num_beamsplitters + inter.num_spaces + inter.num_lenses + inter.num_gratings + inter.num_diodes + inter.num_light_inputs + inter.num_modulators + inter.num_dbss + inter.num_signals; break; case OUT: overall_component_index = inter.num_mirrors + inter.num_beamsplitters + inter.num_spaces + inter.num_lenses + inter.num_gratings + inter.num_diodes + inter.num_light_inputs + inter.num_modulators + inter.num_dbss + inter.num_signals + inter.num_variables; break; case BEAMPARAM: overall_component_index = inter.num_mirrors + inter.num_beamsplitters + inter.num_spaces + inter.num_lenses + inter.num_gratings + inter.num_diodes + inter.num_light_inputs + inter.num_modulators + inter.num_dbss + inter.num_signals + inter.num_variables + inter.num_outputs; break; default: bug_error("no such type"); } return (overall_component_index + component_index); } //! Given the function name return the pointer to the function output /*! * \param function pointer to the function output * \param input_string input string of function name * * \todo can be removed?? * \todo untested * * \see Test_get_function() */ int get_function(double **function, const char *input_string) { *function = NULL; if (input_string == NULL) { return NOT_FOUND; } if (*input_string == '\0') { return NOT_FOUND; } int function_index; for (function_index = 0; function_index < inter.num_func_cmds; function_index++) { if (strcmp(input_string, inter.function_list[function_index].name) == 0) { *function = &inter.function_list[function_index].result; return function_index; } } return NOT_FOUND; } //! Give the index of a set command for a given 'set' name /*! * \param input_string input string of the component to set * * \todo this function not used anywhere in source. Remove??? * * \todo untested * \todo set_nr -> give a better name!!! * * \see Test_set_nr() */ int set_nr(const char *input_string) { if (input_string == NULL) { return NOT_FOUND; } if (*input_string == '\0') { return NOT_FOUND; } int set_cmd_index; for (set_cmd_index = 0; set_cmd_index < inter.num_set_cmds; set_cmd_index++) { if (strcmp(input_string, inter.set_list[set_cmd_index].name) == 0) { return set_cmd_index; } } if (strcasecmp(input_string, "mx3") == 0) { return (MX3); } if (strcasecmp(input_string, "mx2") == 0) { return (MX2); } if (strcasecmp(input_string, "mx1") == 0) { return (MX1); } if (strcasecmp(input_string, "mx") == 0) { return (MX1); } if (strcasecmp(input_string, "x3") == 0) { return (X3); } if (strcasecmp(input_string, "x2") == 0) { return (X2); } if (strcasecmp(input_string, "x1") == 0) { return (X1); } if (strcasecmp(input_string, "x") == 0) { return (X1); } return NOT_FOUND; } //! Given a component name return the component number /*! * \param component_name_string input string of component name * * \see Test_get_component_index_from_name() */ int get_component_index_from_name(const char *component_name_string) { int j; j = 0; if (component_name_string == NULL) { return NOT_FOUND; } if (*component_name_string == '\0') { return NOT_FOUND; } int i; for (i = 0; i < inter.num_mirrors; i++) { if (strcmp(component_name_string, inter.mirror_list[i].name) == 0) { return j + i; } } j += inter.num_mirrors; for (i = 0; i < inter.num_beamsplitters; i++) { if (strcmp(component_name_string, inter.bs_list[i].name) == 0) { return j + i; } } j += inter.num_beamsplitters; for (i = 0; i < inter.num_gratings; i++) { if (strcmp(component_name_string, inter.grating_list[i].name) == 0) { return j + i; } } j += inter.num_gratings; for (i = 0; i < inter.num_spaces; i++) { if (strcmp(component_name_string, inter.space_list[i].name) == 0) { return j + i; } } j += inter.num_spaces; for (i = 0; i < inter.num_sagnacs; i++) { if (strcmp(component_name_string, inter.sagnac_list[i].name) == 0) { return j + i; } } j += inter.num_sagnacs; for (i = 0; i < inter.num_lenses; i++) { if (strcmp(component_name_string, inter.lens_list[i].name) == 0) { return j + i; } } j += inter.num_lenses; for (i = 0; i < inter.num_diodes; i++) { if (strcmp(component_name_string, inter.diode_list[i].name) == 0) { return j + i; } } j += inter.num_diodes; for (i = 0; i < inter.num_light_inputs; i++) { if (strcmp(component_name_string, inter.light_in_list[i].name) == 0) { return j + i; } } j += inter.num_light_inputs; for (i = 0; i < inter.num_modulators; i++) { if (strcmp(component_name_string, inter.modulator_list[i].name) == 0) { return j + i; } } j += inter.num_modulators; for (i = 0; i < inter.num_dbss; i++) { if (strcmp(component_name_string, inter.dbs_list[i].name) == 0) { return j + i; } } j += inter.num_dbss; for (i = 0; i < inter.num_signals; i++) { if (strcmp(component_name_string, inter.signal_list[i].name) == 0) { return j + i; } } j += inter.num_signals; for (i = 0; i < inter.num_variables; i++) { if (strcmp(component_name_string, inter.variable_list[i].name) == 0) { return j + i; } } j += inter.num_variables; for (i = 0; i < inter.num_outputs; i++) { if (strcmp(component_name_string, inter.output_data_list[i].name) == 0) { return j + i; } } j += inter.num_outputs; for (i = 0; i < inter.num_gauss_cmds; i++) { if (strcmp(component_name_string, inter.gauss_list[i].name) == 0) { return j + i; } } j += inter.num_gauss_cmds; for(i=0; i< inter.num_transfer_funcs; i++){ if(strcmp(component_name_string, inter.tf_list[i].name) == 0){ return j+i; } } j += inter.num_transfer_funcs; for(i=0; i< inter.num_locks; i++){ if(strcmp(component_name_string, inter.lock_list[i].name) == 0){ return j+i; } } return NOT_FOUND; } //! Given a component index return the component type /*! * \param component_index the component index * * \see Test_get_component_type() */ int get_component_type(int component_index) { if (component_index == NOT_FOUND) {// experimental: node used only on one side return component_index; } else if (component_index < inter.num_mirrors) { return MIRROR; } else if ((component_index -= inter.num_mirrors) < inter.num_beamsplitters) { return BEAMSPLITTER; } else if ((component_index -= inter.num_beamsplitters) < inter.num_gratings) { return GRATING; } else if ((component_index -= inter.num_gratings) < inter.num_spaces) { return SPACE; } else if ((component_index -= inter.num_spaces) < inter.num_lenses) { return LENS; } else if ((component_index -= inter.num_lenses) < inter.num_diodes) { return DIODE; } else if ((component_index -= inter.num_diodes) < inter.num_light_inputs) { return LIGHT_INPUT; } else if ((component_index -= inter.num_light_inputs) < inter.num_modulators) { return MODULATOR; } else if ((component_index -= inter.num_modulators) < inter.num_dbss) { return DBS; } else if ((component_index -= inter.num_dbss) < inter.num_signals) { return FSIG; } else if ((component_index -= inter.num_signals) < inter.num_variables) { return VARIABLE; } else if ((component_index -= inter.num_variables) < inter.num_outputs) { return OUT; } else if ((component_index -= inter.num_outputs) < inter.num_gauss_cmds) { return BEAMPARAM; } else if ((component_index -= inter.num_gauss_cmds) < inter.num_transfer_funcs) { return TF; } else if ((component_index -= inter.num_transfer_funcs) < inter.num_locks) { return LOCK_CMD; } else { bug_error("unknown type"); } return NOT_FOUND; } //! \brief Given an overall component index return the index with respect //! to only the same type of component /*! * \param overall_component_index the overall component index * * \todo not completely tested * * \see Test_get_component_index() */ int get_component_index(int overall_component_index) { if (overall_component_index == NOT_FOUND) { // experimental: node used only on one side return overall_component_index; } else if (overall_component_index < inter.num_mirrors) { return overall_component_index; } else if ((overall_component_index -= inter.num_mirrors) < inter.num_beamsplitters) { return overall_component_index; } else if ((overall_component_index -= inter.num_beamsplitters) < inter.num_gratings) { return overall_component_index; } else if ((overall_component_index -= inter.num_gratings) < inter.num_spaces) { return overall_component_index; } else if ((overall_component_index -= inter.num_spaces) < inter.num_lenses) { return overall_component_index; } else if ((overall_component_index -= inter.num_lenses) < inter.num_diodes) { return overall_component_index; } else if ((overall_component_index -= inter.num_diodes) < inter.num_light_inputs) { return overall_component_index; } else if ((overall_component_index -= inter.num_light_inputs) < inter.num_modulators) { return overall_component_index; } else if ((overall_component_index -= inter.num_modulators) < inter.num_dbss) { return overall_component_index; } else if ((overall_component_index -= inter.num_dbss) < inter.num_signals) { return overall_component_index; } else if ((overall_component_index -= inter.num_signals) < inter.num_variables) { return overall_component_index; } else if ((overall_component_index -= inter.num_variables) < inter.num_outputs) { return overall_component_index; } else if ((overall_component_index -= inter.num_outputs) < inter.num_gauss_cmds) { return overall_component_index; } else if ((overall_component_index -= inter.num_gauss_cmds) < inter.num_transfer_funcs) { return overall_component_index; } else if ((overall_component_index -= inter.num_transfer_funcs) < inter.num_locks) { return overall_component_index; } else { bug_error("unknown type"); } return NOT_FOUND; } //! \brief Given a component index, return the component type and the index //! is reduced with respect to the type /*! * \param component_index the component index * * \see Test_get_component_type_decriment_index() */ int get_component_type_decriment_index(int *component_index) { if (*component_index >= 0) { if (*component_index < inter.num_mirrors) { return MIRROR; } else if ((*component_index -= inter.num_mirrors) < inter.num_beamsplitters) { return BEAMSPLITTER; } else if ((*component_index -= inter.num_beamsplitters) < inter.num_gratings) { return GRATING; } else if ((*component_index -= inter.num_gratings) < inter.num_spaces) { return SPACE; } else if ((*component_index -= inter.num_spaces) < inter.num_sagnacs) { return SAGNAC; } else if ((*component_index -= inter.num_sagnacs) < inter.num_lenses) { return LENS; } else if ((*component_index -= inter.num_lenses) < inter.num_diodes) { return DIODE; } else if ((*component_index -= inter.num_diodes) < inter.num_light_inputs) { return LIGHT_INPUT; } else if ((*component_index -= inter.num_light_inputs) < inter.num_modulators) { return MODULATOR; } else if ((*component_index -= inter.num_modulators) < inter.num_dbss) { return DBS; } else if ((*component_index -= inter.num_dbss) < inter.num_signals) { return FSIG; } else if ((*component_index -= inter.num_signals) < inter.num_variables) { return VARIABLE; } else if ((*component_index -= inter.num_variables) < inter.num_outputs) { return OUT; } else if ((*component_index -= inter.num_outputs) < inter.num_gauss_cmds) { return BEAMPARAM; } else if ((*component_index -= inter.num_gauss_cmds) < inter.num_transfer_funcs) { return TF; } else if ((*component_index -= inter.num_transfer_funcs) < inter.num_locks) { return LOCK_CMD; } } return NOT_FOUND; } //! Given a component index return the component name as defined in .kat file /*! * \param component_index the component index * * \see Test_get_component_name() */ char *get_component_name(int component_index) { if (component_index < 0) { return duplicate_string("---"); } if (component_index < inter.num_mirrors) { return duplicate_string(inter.mirror_list[component_index].name); } else if ((component_index -= inter.num_mirrors) < inter.num_beamsplitters) { return duplicate_string(inter.bs_list[component_index].name); } else if ((component_index -= inter.num_beamsplitters) < inter.num_gratings) { return duplicate_string(inter.grating_list[component_index].name); } else if ((component_index -= inter.num_gratings) < inter.num_spaces) { return duplicate_string(inter.space_list[component_index].name); } else if ((component_index -= inter.num_spaces) < inter.num_lenses) { return duplicate_string(inter.lens_list[component_index].name); } else if ((component_index -= inter.num_lenses) < inter.num_diodes) { return duplicate_string(inter.diode_list[component_index].name); } else if ((component_index -= inter.num_diodes) < inter.num_light_inputs) { return duplicate_string(inter.light_in_list[component_index].name); } else if ((component_index -= inter.num_light_inputs) < inter.num_modulators) { return duplicate_string(inter.modulator_list[component_index].name); } else if ((component_index -= inter.num_modulators) < inter.num_dbss) { return duplicate_string(inter.dbs_list[component_index].name); } else if ((component_index -= inter.num_dbss) < inter.num_signals) { return duplicate_string(inter.signal_list[component_index].name); } else if ((component_index -= inter.num_signals) < inter.num_variables) { return duplicate_string(inter.variable_list[component_index].name); } else if ((component_index -= inter.num_variables) < inter.num_outputs) { return duplicate_string(inter.output_data_list[component_index].name); } else if ((component_index -= inter.num_outputs) < inter.num_gauss_cmds) { return duplicate_string(inter.gauss_list[component_index].name); } else if ((component_index -= inter.num_gauss_cmds) < inter.num_transfer_funcs) { return duplicate_string(inter.tf_list[component_index].name); } else if ((component_index -= inter.num_transfer_funcs) < inter.num_locks) { return duplicate_string(inter.lock_list[component_index].name); } else { bug_error("component index out of range"); } return duplicate_string("ERROR"); } //! Given a component index return the name of the component type /*! * \param component_index the component index * * \see Test_get_component_type_name() */ char *get_component_type_name(int component_index) { if (component_index < 0) { return duplicate_string("---"); } if (component_index < inter.num_mirrors) { return duplicate_string("MIRROR"); } else if ((component_index -= inter.num_mirrors) < inter.num_beamsplitters) { return duplicate_string("BEAM SPLITTER"); } else if ((component_index -= inter.num_beamsplitters) < inter.num_gratings) { return duplicate_string("GRATING"); } else if ((component_index -= inter.num_gratings) < inter.num_spaces) { return duplicate_string("SPACE"); } else if ((component_index -= inter.num_spaces) < inter.num_lenses) { return duplicate_string("LENS"); } else if ((component_index -= inter.num_lenses) < inter.num_diodes) { return duplicate_string("DIODE"); } else if ((component_index -= inter.num_diodes) < inter.num_light_inputs) { return duplicate_string("LIGHT INPUT"); } else if ((component_index -= inter.num_light_inputs) < inter.num_modulators) { return duplicate_string("MODULATOR"); } else if ((component_index -= inter.num_modulators) < inter.num_dbss) { return duplicate_string("DBS"); } else if ((component_index -= inter.num_dbss) < inter.num_signals) { return duplicate_string("SIGNAL"); } else if ((component_index -= inter.num_signals) < inter.num_variables) { return duplicate_string("VARIABLE"); } else if ((component_index -= inter.num_variables) < inter.num_outputs) { return duplicate_string("OUTPUT"); } else if ((component_index -= inter.num_outputs) < inter.num_gauss_cmds) { return duplicate_string("BEAM PARAMETER"); } else if ((component_index -= inter.num_gauss_cmds) < inter.num_transfer_funcs) { return duplicate_string("TF"); } else if ((component_index -= inter.num_transfer_funcs) < inter.num_locks) { return duplicate_string("LOCK"); } else { bug_error("component index out of range"); } return duplicate_string("ERROR"); } //! Given a detector name return the detector number /*! * \param detector_name input string of the detector name * * \see Test_get_detector_index_from_name() */ int get_detector_or_node_index_from_name(const char *detector_name) { if (detector_name == NULL) { return NOT_FOUND; } if (*detector_name == '\0') { return NOT_FOUND; } int j = -1; int i; for (i = 0; i < inter.num_outputs; i++) { if (strcmp(detector_name, inter.output_data_list[i].name) == 0) { return (j - i); } } j = j - inter.num_outputs; for (i = 0; i < inter.num_gauss_cmds; i++) { if (strcmp(detector_name, inter.gauss_list[i].name) == 0) { return (j - i); } } return NOT_FOUND; } //! Given a node name return the node number /*! * \param node_name input string of the node name * */ int get_node_index_from_name(const char *node_name) { if (node_name == NULL) { return NOT_FOUND; } if (*node_name == '\0') { return NOT_FOUND; } int node_index; for (node_index = 0; node_index < inter.num_nodes; node_index++) { if (strcmp(node_name, inter.node_list[node_index].name) == 0) { return (node_index); } } return NOT_FOUND; } //! Given a detector number return the detector name /*! * \param detector_index detector index * * \todo function not used in code: remove??? * * \see Test_get_detector_name() */ char *get_detector_name(int detector_index) { if (detector_index > 0) { bug_error("detector index greater than zero"); } else { detector_index = -detector_index - 1; } if (detector_index < inter.num_outputs) { return duplicate_string(inter.output_data_list[detector_index].name); } else if (detector_index -= inter.num_outputs < inter.num_gauss_cmds) { return duplicate_string(inter.gauss_list[detector_index].name); } else { bug_error("detector/node index out of range"); } return duplicate_string("ERROR"); } //! Convert character to integer /*! * \param c the character to convert * \param i the integer value to return * * \todo this function isn't used in the code: remove??? * * \todo this should be updated to a better implementation * * \see Test_ctoi() */ int ctoi(char c, int *i) { if (c < 48 || c > 57) { return (1); } else { *i = (int) c - 48; } return (0); } //! My version of ascii to integer /*! * \param input_string ascii string to convert * \param converted_integer pointer to returned integer * * \todo maybe rename to something more intuitive * \see Test_my_atoi() */ int my_atoi(char *input_string, int *converted_integer) { char xxx[LINE_LEN] = {0}; char *nptr; nptr = xxx; *converted_integer = strtol(input_string, &nptr, 10); if (nptr == NULL) { return (0); } else { return (strlen(nptr)); } } //! Convert ascii characters to double precision numbers /*! * \param input_string ascii string to convert * \param converted_double pointer to returned double precision number * * \return true for success, false for failure * * \see Test_atod() */ bool atod(const char *input_string, double *converted_double) { char *magnitude_suffix_string; *converted_double = strtod(input_string, &magnitude_suffix_string); if (magnitude_suffix_string == input_string) { /* error */ if (strncmp(input_string, "Inf", 3)) { // switched of warning because this makes the parsing more complex and // provides not much extra info, except for debugging if (inter.debug & 512) { message("atod: couldn't read number from '%s'\n", input_string); } return true; } else { /* sollte das untere sein aber kann warnung nicht mehr sehen! */ *converted_double = 0.0; /* *d=1./0.; */ *magnitude_suffix_string = '\0'; } } if (*magnitude_suffix_string == 'f') { *converted_double *= 1e-15; ++magnitude_suffix_string; } else if (*magnitude_suffix_string == 'p') { *converted_double *= 1e-12; ++magnitude_suffix_string; } else if (*magnitude_suffix_string == 'n') { *converted_double *= 1e-9; ++magnitude_suffix_string; } else if (*magnitude_suffix_string == 'u') { *converted_double *= 1e-6; ++magnitude_suffix_string; } else if (*magnitude_suffix_string == 'm') { *converted_double *= 1e-3; ++magnitude_suffix_string; } else if (*magnitude_suffix_string == 'k') { *converted_double *= 1e3; ++magnitude_suffix_string; } else if (*magnitude_suffix_string == 'M') { *converted_double *= 1e6; ++magnitude_suffix_string; } else if (*magnitude_suffix_string == 'G') { *converted_double *= 1e9; ++magnitude_suffix_string; } /* Now the rest should be only whitespace */ while (*magnitude_suffix_string == ' ' || *magnitude_suffix_string == '\t' || *magnitude_suffix_string == '\n') { ++magnitude_suffix_string; /* skip whitespace */ } if (*magnitude_suffix_string != '\0' && *magnitude_suffix_string != '#') { warn("In number '%s': Don't know how to handle '%s'\n", input_string, magnitude_suffix_string); } return false; } //! Construct output form for printing complex numbers /*! * \param z complex number to output * * \see Test_complex_form() */ char *complex_form(const complex_t z) { static char result[FORM_LEN] = {0}; if (fabs(z.re) < init.small && fabs(z.im) < init.small) { sprintf(result, "0"); } else if (fabs(z.re) >= init.small && fabs(z.im) < init.small) { sprintf(result, "%g", z.re); } else if (fabs(z.re) < init.small && fabs(z.im) >= init.small) { sprintf(result, "%gj", z.im); } else { if (z.im < 0) { sprintf(result, "%g-%gj", z.re, -z.im); } else { sprintf(result, "%g+%gj", z.re, z.im); } } return duplicate_string(result); } //! \brief Construct output form for printing complex numbers but with //! 15 decimal places of precision /*! * \param z complex number to output * * \see Test_complex_form15() */ char *complex_form15(const complex_t z) { static char result[FORM_LEN] = {0}; if (fabs(z.re) < init.small && fabs(z.im) < init.small) { sprintf(result, "0"); } else if (fabs(z.re) >= init.small && fabs(z.im) < init.small) { sprintf(result, "%.15g + 0i", z.re); } else if (fabs(z.re) < init.small && fabs(z.im) >= init.small) { if (z.im < 0) { sprintf(result, "0 - %.15gi", -z.im); } else { sprintf(result, "0 + %.15gi", z.im); } } else { if (z.im < 0) { sprintf(result, "(%.15g - %.15gi)", z.re, -z.im); } else { sprintf(result, "(%.15g + %.15gi)", z.re, z.im); } } return duplicate_string(result); } //! Frequency output format /*! * \param f frequency to output * * \see Test_freq_form() */ char *freq_form(double f) { static char result[FORM_LEN] = {0}; char prefix[2] = {0}; si_prefix(&f, prefix); sprintf(result, "%.15g %sHz", f, prefix); return duplicate_string(result); } //! Integer output format /*! * \param output_integer integer number to output * * \see Test_int_form() */ char *int_form(int output_integer) { static char result[FORM_LEN] = {0}; sprintf(result, "%d", output_integer); return duplicate_string(result); } //! Double output using an SI prefix /*! * \param x input double precision number to format * * \see Test_double_form() */ char *double_form(double x) { static char result[FORM_LEN] = {0}; char prefix[2] = {0}; double xa = fabs(x); int sign = (x < 0) ? -1 : 1; si_prefix(&xa, prefix); sprintf(result, "%.4g%s", xa * sign, prefix); return duplicate_string(result); } //! Double output using an SI prefix; more precision than double_form /*! * In other words, extended double form * * \see double_form * * \param x Input double precision number to format * * \see Test_xdouble_form() */ char *xdouble_form(double x) { /* higher precision */ static char result[FORM_LEN] = {0}; char prefix[2] = {0}; double xa = fabs(x); int sign = (x < 0) ? -1 : 1; si_prefix(&xa, prefix); if(options.perl1) sprintf(result, "%.15g%s", xa * sign, prefix); else sprintf(result, "%.8g%s", xa * sign, prefix); return duplicate_string(result); } //! Construct and return the long output format string /*! * \param value_to_print the number to be printed * \param padding_amount the amount of preceding space in the number * * \see Test_long_form() */ char *long_form(double value_to_print, int padding_amount) { static char result[FORM_LEN] = {0}; if (fabs(value_to_print) < SMALL) { value_to_print = 0; } switch (padding_amount) { case 1: sprintf(result, "%1g ", value_to_print); break; case 2: sprintf(result, "%2g ", value_to_print); break; case 3: sprintf(result, "%3g ", value_to_print); break; case 4: sprintf(result, "%4g ", value_to_print); break; case 5: sprintf(result, "%5g ", value_to_print); break; case 6: sprintf(result, "%6g ", value_to_print); break; case 7: sprintf(result, "%7g ", value_to_print); break; case 8: sprintf(result, "%8g ", value_to_print); break; case 9: sprintf(result, "%9g ", value_to_print); break; case 10: sprintf(result, "%10g ", value_to_print); break; case 11: sprintf(result, "%11g ", value_to_print); break; case 12: sprintf(result, "%12g ", value_to_print); break; default: sprintf(result, "%g ", value_to_print); } return duplicate_string(result); } //! Generate the SI prefix for numbers /*! * \param double_to_interpret the number for which to determine the SI prefix * \param prefix_string the SI prefix to use (returned) * * \see Test_si_prefix() */ void si_prefix(double *double_to_interpret, char *prefix_string) { if (*double_to_interpret == 0) { strcpy(prefix_string, ""); return; } else if (fabs(*double_to_interpret) < 1e-12) { *double_to_interpret *= 1e15; strcpy(prefix_string, "f"); return; } else if (fabs(*double_to_interpret) < 1e-9) { *double_to_interpret *= 1e12; strcpy(prefix_string, "p"); return; } else if (fabs(*double_to_interpret) < 1e-6) { *double_to_interpret *= 1e9; strcpy(prefix_string, "n"); return; } else if (fabs(*double_to_interpret) < 1e-3) { *double_to_interpret *= 1e6; strcpy(prefix_string, "u"); return; } else if (fabs(*double_to_interpret) < 1) { *double_to_interpret *= 1e3; strcpy(prefix_string, "m"); return; } else if (fabs(*double_to_interpret) < 1e3) { strcpy(prefix_string, ""); return; } else if (fabs(*double_to_interpret) < 1e6) { *double_to_interpret *= 1e-3; strcpy(prefix_string, "k"); return; } else if (fabs(*double_to_interpret) < 1e9) { *double_to_interpret *= 1e-6; strcpy(prefix_string, "M"); return; } *double_to_interpret *= 1e-9; strcpy(prefix_string, "G"); return; } // Bessel functions from Takuya Ooura // (see http://momonga.t.u-tokyo.ac.jp/~ooura/). Compared to // available tables the old version differed in the 7th/8th // significant digit whereas the new is correct to the accuracy of // the table in `Handbook of Mathematical Functions...' // by Abromowitz, Stegun (J0: 15 digits, J1, J2: 10 digits). //! Zeroth order Bessel function of the first kind /*! * \param x position on x-axis * * \see Test_dbesj0() */ double dbesj0(double x) { int k; double w, t, y, v, theta; static double a[8] = { -2.3655394e-12, 4.70889868e-10, -6.78167892231e-8, 6.7816840038636e-6, -4.340277777716935e-4, 0.0156249999999992397, -0.2499999999999999638, 0.9999999999999999997 }; static double b[65] = { 6.26681117e-11, -2.2270614428e-9, 6.62981656302e-8, -1.6268486502196e-6, 3.21978384111685e-5, -5.00523773331583e-4, 0.0059060313537449816, -0.0505265323740109701, 0.2936432097610503985, -1.0482565081091638637, 1.9181123286040428113, -1.13191994752217001, -0.1965480952704682, 4.57457332e-11, -1.5814772025e-9, 4.55487446311e-8, -1.0735201286233e-6, 2.02015179970014e-5, -2.942392368203808e-4, 0.0031801987726150648, -0.0239875209742846362, 0.1141447698973777641, -0.2766726722823530233, 0.1088620480970941648, 0.5136514645381999197, -0.2100594022073706033, 3.31366618e-11, -1.1119090229e-9, 3.08823040363e-8, -6.956602653104e-7, 1.23499947481762e-5, -1.66295194539618e-4, 0.0016048663165678412, -0.0100785479932760966, 0.0328996815223415274, -0.0056168761733860688, -0.2341096400274429386, 0.2551729256776404262, 0.2288438186148935667, 2.38007203e-11, -7.731046439e-10, 2.06237001152e-8, -4.412291442285e-7, 7.3107766249655e-6, -8.91749801028666e-5, 7.34165451384135e-4, -0.0033303085445352071, 0.0015425853045205717, 0.0521100583113136379, -0.1334447768979217815, -0.1401330292364750968, 0.2685616168804818919, 1.6935595e-11, -5.308092192e-10, 1.35323005576e-8, -2.726650587978e-7, 4.151324014176e-6, -4.43353052220157e-5, 2.815740758993879e-4, -4.393235121629007e-4, -0.0067573531105799347, 0.0369141914660130814, 0.0081673361942996237, -0.257338128589888186, 0.0459580257102978932 }; static double c[70] = { -3.009451757e-11, -1.4958003844e-10, 5.06854544776e-9, 1.863564222012e-8, -6.0304249068078e-7, -1.47686259937403e-6, 4.714331342682714e-5, 6.286305481740818e-5, -0.00214137170594124344, -8.9157336676889788e-4, 0.04508258728666024989, -0.00490362805828762224, -0.27312196367405374426, 0.04193925184293450356, -7.1245356e-12, -4.1170814825e-10, 1.38012624364e-9, 5.704447670683e-8, -1.9026363528842e-7, -5.33925032409729e-6, 1.736064885538091e-5, 3.0692619152608375e-4, -9.2598938200644367e-4, -0.00917934265960017663, 0.02287952522866389076, 0.10545197546252853195, -0.16126443075752985095, -0.19392874768742235538, 2.128344556e-11, -3.1053910272e-10, -3.34979293158e-9, 4.50723289505e-8, 3.6437959146427e-7, -4.46421436266678e-6, -2.523429344576552e-5, 2.7519882931758163e-4, 9.7185076358599358e-4, -0.00898326746345390692, -0.01665959196063987584, 0.11456933464891967814, 0.07885001422733148815, -0.23664819446234712621, 3.035295055e-11, 5.486066835e-11, -5.01026824811e-9, -5.0124684786e-9, 5.8012340163034e-7, 1.6788922416169e-7, -4.373270270147275e-5, 1.183898532719802e-5, 0.00189863342862291449, -0.0011375924956163613, -0.03846797195329871681, 0.02389746880951420335, 0.22837862066532347461, -0.06765394811166522844, 1.279875977e-11, 3.5925958103e-10, -2.28037105967e-9, -4.852770517176e-8, 2.8696428000189e-7, 4.40131125178642e-6, -2.366617753349105e-5, -2.4412456252884129e-4, 0.00113028178539430542, 0.0070847051391978908, -0.02526914792327618386, -0.08006137953480093426, 0.16548380461475971846, 0.14688405470042110229 }; static double d[52] = { 1.059601355592185731e-14, -2.71150591218550377e-13, 8.6514809056201638e-12, -4.6264028554286627e-10, 5.0815403835647104e-8, -1.76722552048141208e-5, 0.16286750396763997378, 2.949651820598278873e-13, -8.818215611676125741e-12, 3.571119876162253451e-10, -2.63192412099371706e-8, 4.709502795656698909e-6, -0.005208333333333283282, 7.18344107717531977e-15, -2.51623725588410308e-13, 8.6017784918920604e-12, -4.6256876614290359e-10, 5.0815343220437937e-8, -1.7672255176494197e-5, 0.16286750396763433767, 2.2327570859680094777e-13, -8.464594853517051292e-12, 3.563766464349055183e-10, -2.631843986737892965e-8, 4.70950234228865941e-6, -0.0052083333332278466225, 5.15413392842889366e-15, -2.27740238380640162e-13, 8.4827767197609014e-12, -4.6224753682737618e-10, 5.0814848128929134e-8, -1.7672254763876748e-5, 0.16286750396748926663, 1.7316195320192170887e-13, -7.971122772293919646e-12, 3.544039469911895749e-10, -2.631443902081701081e-8, 4.709498228695400603e-6, -0.005208333331514365361, 3.84653681453798517e-15, -2.04464520778789011e-13, 8.3089298605177838e-12, -4.6155016158412096e-10, 5.081326369646665e-8, -1.76722528311426167e-5, 0.1628675039665006593, 1.3797879972460878797e-13, -7.448089381011684812e-12, 3.51273379710695978e-10, -2.630500895563592722e-8, 4.709483934775839193e-6, -0.0052083333227940760113 }; w = fabs(x); if (w < 1) { t = w * w; y = ((((((a[0] * t + a[1]) * t + a[2]) * t + a[3]) * t + a[4]) * t + a[5]) * t + a[6]) * t + a[7]; } else if (w < 8.5) { t = w * w * 0.0625; k = (int) t; t -= k + 0.5; k *= 13; y = (((((((((((b[k] * t + b[k + 1]) * t + b[k + 2]) * t + b[k + 3]) * t + b[k + 4]) * t + b[k + 5]) * t + b[k + 6]) * t + b[k + 7]) * t + b[k + 8]) * t + b[k + 9]) * t + b[k + 10]) * t + b[k + 11]) * t + b[k + 12]; } else if (w < 12.5) { k = (int) w; t = w - (k + 0.5); k = 14 * (k - 8); y = ((((((((((((c[k] * t + c[k + 1]) * t + c[k + 2]) * t + c[k + 3]) * t + c[k + 4]) * t + c[k + 5]) * t + c[k + 6]) * t + c[k + 7]) * t + c[k + 8]) * t + c[k + 9]) * t + c[k + 10]) * t + c[k + 11]) * t + c[k + 12]) * t + c[k + 13]; } else { v = 24 / w; t = v * v; k = 13 * ((int) t); y = ((((((d[k] * t + d[k + 1]) * t + d[k + 2]) * t + d[k + 3]) * t + d[k + 4]) * t + d[k + 5]) * t + d[k + 6]) * sqrt(v); theta = (((((d[k + 7] * t + d[k + 8]) * t + d[k + 9]) * t + d[k + 10]) * t + d[k + 11]) * t + d[k + 12]) * v - 0.78539816339744830962; y *= cos(w + theta); } return y; } //! First order Bessel function of the first kind /*! * \param x position on x-axis * * \todo else branch is untested * * \see Test_dbesj1() */ double dbesj1(double x) { int k; double w, t, y, v, theta; static double a[8] = { -1.4810349e-13, 3.363594618e-11, -5.65140051697e-9, 6.7816840144764e-7, -5.425347222188379e-5, 0.00260416666666662438, -0.06249999999999999799, 0.49999999999999999998 }; static double b[65] = { 2.43721316e-12, -9.400554763e-11, 3.0605338998e-9, -8.287270492518e-8, 1.83020515991344e-6, -3.219783841164382e-5, 4.3795830161515318e-4, -0.00442952351530868999, 0.03157908273375945955, -0.14682160488052520107, 0.39309619054093640008, -0.4795280821510107028, 0.1414899934402712514, 1.82119257e-12, -6.862117678e-11, 2.1732790836e-9, -5.69359291782e-8, 1.20771046483277e-6, -2.020151799736374e-5, 2.5745933218048448e-4, -0.00238514907946126334, 0.01499220060892984289, -0.05707238494868888345, 0.10375225210588234727, -0.02721551202427354117, -0.06420643306727498985, 1.352611196e-12, -4.9706947875e-11, 1.527944986332e-9, -3.8602878823401e-8, 7.82618036237845e-7, -1.23499947484511e-5, 1.45508295194426686e-4, -0.001203649737425854162, 0.006299092495799005109, -0.016449840761170764763, 0.002106328565019748701, 0.05852741000686073465, -0.031896615709705053191, 9.97982124e-13, -3.5702556073e-11, 1.062332772617e-9, -2.5779624221725e-8, 4.96382962683556e-7, -7.310776625173004e-6, 7.8028107569541842e-5, -5.50624088538081113e-4, 0.002081442840335570371, -7.71292652260286633e-4, -0.019541271866742634199, 0.033361194224480445382, 0.017516628654559387164, 7.31050661e-13, -2.5404499912e-11, 7.29360079088e-10, -1.6915375004937e-8, 3.06748319652546e-7, -4.151324014331739e-6, 3.8793392054271497e-5, -2.11180556924525773e-4, 2.74577195102593786e-4, 0.003378676555289966782, -0.013842821799754920148, -0.002041834048574905921, 0.032167266073736023299 }; static double c[70] = { -1.185964494e-11, 3.9110295657e-10, 1.80385519493e-9, -5.575391345723e-8, -1.8635897017174e-7, 5.42738239401869e-6, 1.181490114244279e-5, -3.300031939852107e-4, -3.7717832892725053e-4, 0.01070685852970608288, 0.00356629346707622489, -0.13524776185998074716, 0.00980725611657523952, 0.27312196367405374425, -3.029591097e-11, 9.259293559e-11, 4.96321971223e-9, -1.518137078639e-8, -5.7045127595547e-7, 1.71237271302072e-6, 4.271400348035384e-5, -1.2152454198713258e-4, -0.00184155714921474963, 0.00462994691003219055, 0.03671737063840232452, -0.06863857568599167175, -0.21090395092505707655, 0.16126443075752985095, -2.19760208e-11, -2.7659100729e-10, 3.74295124827e-9, 3.684765777023e-8, -4.5072801091574e-7, -3.27941630669276e-6, 3.5713715545163e-5, 1.7664005411843533e-4, -0.00165119297594774104, -0.00485925381792986774, 0.03593306985381680131, 0.04997877588191962563, -0.22913866929783936544, -0.07885001422733148814, 5.16292316e-12, -3.9445956763e-10, -6.6220021263e-10, 5.511286218639e-8, 5.01257940078e-8, -5.22111059203425e-6, -1.34311394455105e-6, 3.0612891890766805e-4, -7.103391195326182e-5, -0.00949316714311443491, 0.00455036998246516948, 0.11540391585989614784, -0.04779493761902840455, -0.2283786206653234746, 2.697817493e-11, -1.6633326949e-10, -4.3313486035e-9, 2.508404686362e-8, 4.8528284780984e-7, -2.58267851112118e-6, -3.521049080466759e-5, 1.6566324273339952e-4, 0.00146474737522491617, -0.00565140892697147306, -0.028338820556793004, 0.07580744376982855057, 0.16012275906960187978, -0.16548380461475971845 }; static double d[52] = { -1.272346002224188092e-14, 3.370464692346669075e-13, -1.144940314335484869e-11, 6.863141561083429745e-10, -9.491933932960924159e-8, 5.301676561445687562e-5, 0.162867503967639974, -3.652982212914147794e-13, 1.151126750560028914e-11, -5.165585095674343486e-10, 4.657991250060549892e-8, -1.186794704692706504e-5, 0.01562499999999994026, -8.713069680903981555e-15, 3.140780373478474935e-13, -1.139089186076256597e-11, 6.862299023338785566e-10, -9.491926788274594674e-8, 5.301676558106268323e-5, 0.162867503967646622, -2.792555727162752006e-13, 1.108650207651756807e-11, -5.156745588549830981e-10, 4.657894859077370979e-8, -1.186794650130550256e-5, 0.01562499999987299901, -6.304859171204770696e-15, 2.857249044208791652e-13, -1.124956921556753188e-11, 6.858482894906716661e-10, -9.49186795351689846e-8, 5.301676509057781574e-5, 0.1628675039678191167, -2.185193490132496053e-13, 1.048820673697426074e-11, -5.132819367467680132e-10, 4.65740943737299422e-8, -1.186794150862988921e-5, 0.01562499999779270706, -4.74041720979200985e-15, 2.578715253644144182e-13, -1.104148898414138857e-11, 6.850134201626289183e-10, -9.49167823417491964e-8, 5.301676277588728159e-5, 0.1628675039690033136, -1.75512205749384229e-13, 9.848723331445182397e-12, -5.094535425482245697e-10, 4.656255982268609304e-8, -1.186792402114394891e-5, 0.01562499998712198636 }; w = fabs(x); if (w < 1) { t = w * w; y = (((((((a[0] * t + a[1]) * t + a[2]) * t + a[3]) * t + a[4]) * t + a[5]) * t + a[6]) * t + a[7]) * w; } else if (w < 8.5) { t = w * w * 0.0625; k = (int) t; t -= k + 0.5; k *= 13; y = ((((((((((((b[k] * t + b[k + 1]) * t + b[k + 2]) * t + b[k + 3]) * t + b[k + 4]) * t + b[k + 5]) * t + b[k + 6]) * t + b[k + 7]) * t + b[k + 8]) * t + b[k + 9]) * t + b[k + 10]) * t + b[k + 11]) * t + b[k + 12]) * w; } else if (w < 12.5) { k = (int) w; t = w - (k + 0.5); k = 14 * (k - 8); y = ((((((((((((c[k] * t + c[k + 1]) * t + c[k + 2]) * t + c[k + 3]) * t + c[k + 4]) * t + c[k + 5]) * t + c[k + 6]) * t + c[k + 7]) * t + c[k + 8]) * t + c[k + 9]) * t + c[k + 10]) * t + c[k + 11]) * t + c[k + 12]) * t + c[k + 13]; } else { v = 24 / w; t = v * v; k = 13 * ((int) t); y = ((((((d[k] * t + d[k + 1]) * t + d[k + 2]) * t + d[k + 3]) * t + d[k + 4]) * t + d[k + 5]) * t + d[k + 6]) * sqrt(v); theta = (((((d[k + 7] * t + d[k + 8]) * t + d[k + 9]) * t + d[k + 10]) * t + d[k + 11]) * t + d[k + 12]) * v - 0.78539816339744830962; y *= sin(w + theta); } return x < 0 ? -y : y; } //! Returns z to the power of x /*! * returns \f$z^x\f$ by using \f$z^x = \exp(\log(z^x)) = \exp(x\log(z))\f$ * * \param z the complex base to be raised * \param x the power to raise the complex number by * * \todo untested * * \see Test_my_cpow() */ complex_t my_cpow(complex_t z, double x) { complex_t z1; z1 = my_clog(z); z1.re *= x; z1.im *= x; return (my_cexp(z1)); } //! Returns the exponent of the complex number z /*! * returns \f$\exp(z)\f$ as * \f$\exp(z) = \exp(x + iy) = \exp(x)(\cos(y) + i\sin(y))\f$ * * \param z input complex number * * \see Test_my_cexp() */ complex_t my_cexp(complex_t z) { double x; complex_t z1; x = exp(z.re); z1.re = x * cos(z.im); z1.im = x * sin(z.im); return (z1); } //! Returns logarithm of the complex number z /*! * returns \f$\log(z)\f$ as * \f$\log(z) = \log(r \exp(i\phi)) = \log(r) + i\phi\f$ * * \param z input complex number * * \see Test_my_clog() */ complex_t my_clog(complex_t z) { double x; complex_t z1; x = sqrt(z.re * z.re + z.im * z.im); z1.re = log(x); z1.im = atan2(z.im, z.re); return (z1); } /** * returns i^k */ complex_t pow_complex_i(int k){ switch(k % 4){ case 0: return complex_1; case 1: case -3: return complex_i; case 2: case -2: return complex_m1; case 3: case -1: return complex_mi; default: bug_error("value not expected"); } return complex_0; } inline int pow_two(int n){ return 1 << n; } inline int pow_neg_1(int n){ return (GSL_IS_EVEN(n) ? 1.0 : -1.0); } //! Calculates in general z^(nom/denom) /*! * \param z input complex number * \param nom nominator of exponent * \param denom denominator of exponent * * \see Test_pow_complex() */ complex_t pow_complex(complex_t z, int nom, int denom) { complex_t z1; double x; if (denom == 0) { bug_error("denominator is zero"); } if (nom == 0) { return (complex_1); } x = (double) nom / (double) denom; if (ceq(z, complex_0)) { if (x > 0) { return (complex_0); } else { bug_error("division by zero (1)"); } } if (nom == denom) { // nom=denom -> result=z return (z); } if (eq(x, vlocal.x_p1) && ceq(z, vlocal.z_in_p1)) { return (vlocal.z_out_p1); } if (eq(x, vlocal.x_p2) && ceq(z, vlocal.z_in_p2)) { return (vlocal.z_out_p2); } vlocal.x_p2 = vlocal.x_p1; vlocal.x_p1 = x; vlocal.z_in_p2 = vlocal.z_in_p1; vlocal.z_in_p1 = z; vlocal.z_out_p2 = vlocal.z_out_p1; z1 = my_clog(z); z1.re *= x; z1.im *= x; vlocal.z_out_p1 = my_cexp(z1); return (vlocal.z_out_p1); } /* ************************************************************************ * C math library * function ZEROIN - obtain a function zero within the given range * * Input * double zeroin(ax, bx, f, tol) * double ax; Root will be seeked for within * double bx; a range [ax, bx] * double (*f)(double x); Name of the function whose zero * will be seeked for * double tol; Acceptable tolerance for the root * value. * May be specified as 0.0 to cause * the program to find the root as * accurate as possible * * Output * Zeroin returns an estimate for the root with accuracy * 4*EPSILON*abs(x) + tol * * Algorithm * G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical * computations. M., Mir, 1980, p.180 of the Russian edition * * The function makes use of the bissection procedure combined with * the linear or quadric inverse interpolation. * At every step program operates on three abscissae - a, b, and c. * b - the last and the best approximation to the root * a - the last but one approximation * c - the last but one or even earlier approximation than a that * 1) |f(b)| <= |f(c)| * 2) f(b) and f(c) have opposite signs, i.e. b and c confine * the root * At every step Zeroin selects one of the two new approximations, the * former being obtained by the bissection procedure and the latter * resulting in the interpolation (if a, b, and c are all different * the quadric interpolation is utilized, otherwise the linear one). * If the latter (i.e. obtained by the interpolation) point is * reasonable (i.e. lies within the current interval [b, c] not being * too close to the boundaries) it is accepted. The bissection result * is used in the other case. Therefore, the range of uncertainty is * ensured to be reduced at least by the factor 1.6 * ************************************************************************ */ #define ZEROEPSILON DBL_EPSILON //!< Roughly close to zero??? //! Obtain a function zero within the given range /*! * \param ax left border of the range in which the root is sought * \param bx right border of the range in which the root is sought * \param f function under investigation * \param tol acceptable tolerance * \return an estimate to the root * * \todo untested * * \see Test_zeroin() */ double zeroin(double ax, double bx, double (*f) (double), double tol) /* Output : An estimate to the root */ /* ax : Left border | of the range */ /* bx : Right border| the root is seeked */ /* f : Function under investigation */ /* tol : Acceptable tolerance */ { double a, b, c; /* Abscissae, descr. see above */ double fa; /* f(a) */ double fb; /* f(b) */ double fc; /* f(c) */ a = ax; b = bx; fa = (*f) (a); fb = (*f) (b); c = a; fc = fa; for (;;) { /* Main iteration loop */ double prev_step = b - a; /* Distance from the last but one */ /* to the last approximation */ double tol_act; /* Actual tolerance */ double p; /* Interpolation step is calcu- */ double q; /* lated in the form p/q; divi- */ /* sion operations is delayed */ /* until the last moment */ double new_step; /* Step at this iteration */ if (fabs(fc) < fabs(fb)) { /* Swap data for b to be the */ a = b; b = c; c = a; /* best approximation */ fa = fb; fb = fc; fc = fa; } tol_act = 2 * ZEROEPSILON * fabs(b) + tol / 2; new_step = (c - b) / 2; if (fabs(new_step) <= tol_act || fb == (double) 0) { return b; /* Acceptable approx. is found */ } /* Decide if the interpolation can be tried */ if (fabs(prev_step) >= tol_act /* If prev_step was large enough */ && fabs(fa) > fabs(fb)) { /* and was in true direction, */ /* Interpolatiom may be tried */ register double t1, cb, t2; /* If we have only two distinct points linear interpolation * can only be applied */ cb = c - b; if (a == c) { t1 = fb / fa; p = cb * t1; q = 1.0 - t1; } else { /* Quadric inverse interpolation */ q = fa / fc; t1 = fb / fc; t2 = fb / fa; p = t2 * (cb * q * (q - t1) - (b - a) * (t1 - 1.0)); q = (q - 1.0) * (t1 - 1.0) * (t2 - 1.0); } if (p > (double) 0) { /* p was calculated with the op- */ q = -q; /* posite sign; make p positive */ } else { /* and assign possible minus to q */ p = -p; } /* If b+p/q falls in [b, c] */ if (p < (0.75 * cb * q - fabs(tol_act * q) / 2) && p < fabs(prev_step * q / 2)) { /* and isn't too large */ new_step = p / q; /* it is accepted */ } /* If p/q is too large then the */ /* bissection procedure can */ /* reduce [b, c] range to more */ /* extent */ } /* Adjust the step to be not less than tolerance */ if (fabs(new_step) < tol_act) { if (new_step > (double) 0) { new_step = tol_act; } else { new_step = -tol_act; } } a = b; fa = fb; /* Save the previous approx. */ b += new_step; fb = (*f) (b); /* Do step to a new approxim. */ /* Adjust c for it to have a sign */ if ((fb > 0 && fc > 0) || (fb < 0 && fc < 0)) { c = a; fc = fa; /* opposite to that of b */ } } } //! Does the "match string" match any part of the input string? /*! * \param input_string The input string to test * \param match_string The string to match within the input string */ bool string_matches(char *input_string, char *match_string) { if (strncasecmp(input_string, match_string, strlen(match_string)) == 0) { return true; } else { return false; } } //! Does the "match string" matches exactly the input string? /*! * \param input_string The input string to test * \param match_string The string to match within the input string */ bool string_matches_exactly(const char *input_string, const char *match_string) { if (strcasecmp(input_string, match_string) == 0) { return true; } else { return false; } } //! Replace substring in string with other string char *replace_str(char *str, char *orig, char *rep) { static char buffer[4096]; char *p; if(!(p = strstr(str, orig))) // Is substring present in string? return str; strncpy(buffer, str, p-str); buffer[p-str] = '\0'; sprintf(buffer+(p-str), "%s%s", rep, p+strlen(orig)); return buffer; } int getNumCores() { #if defined(OSWIN) && !defined(__CYGWIN__) SYSTEM_INFO sysinfo; GetSystemInfo(&sysinfo); return sysinfo.dwNumberOfProcessors; #elif OS==MACOS int nm[2]; size_t len = 4; uint32_t count; nm[0] = CTL_HW; nm[1] = HW_AVAILCPU; sysctl(nm, 2, &count, &len, NULL, 0); if (count < 1) { nm[1] = HW_NCPU; sysctl(nm, 2, &count, &len, NULL, 0); if (count < 1) { count = 1; } } return count; #else return sysconf(_SC_NPROCESSORS_ONLN); #endif } /** * Check to see if Cygwin server is running. If not then the shared memory commands * that the integrators use will not work. * * Other platforms always default to return true * * @return true or false depending if process is running or not */ bool isCygserverRunning() { // check if we have built with SHMGET, if we use shmget then checking if // cygserver is present is necessary #if __CYGWIN__ && HAVE_SHMGET #pragma message "Cuba has been built using shmget(). This will not run on windows machines without cygserver running. See http://kvasir.sr.bham.ac.uk/redmine/issues/111" PROCESSENTRY32 entry; entry.dwSize = sizeof(PROCESSENTRY32); HANDLE snapshot = CreateToolhelp32Snapshot(TH32CS_SNAPPROCESS, (DWORD)NULL); if (Process32First(snapshot, &entry) == TRUE) { while (Process32Next(snapshot, &entry) == TRUE) { if (stricmp(entry.szExeFile, "cygserver.exe") == 0) { HANDLE hProcess = OpenProcess(PROCESS_ALL_ACCESS, FALSE, entry.th32ProcessID); CloseHandle(hProcess); CloseHandle(snapshot); return true; } } } CloseHandle(snapshot); return false; #else return true; // Always return true on other platforms as they can run the various shmget commands #endif } /** * Computes an MD5 hash of some file * * @param mapfile File to compute hash of, should be opened * @param hash Previously allocated array to store hash in, */ void file_MD5(FILE *mapfile, unsigned char *hash) { assert(mapfile != NULL); assert(hash != NULL); MD5_CTX ctx; MD5Init(&ctx); unsigned int filesize; unsigned char *buffer; //calculate file size fseek(mapfile,0, SEEK_END); filesize = ftell(mapfile); fseek(mapfile,0, SEEK_SET); //allocate new memory for storing entire file buffer = (unsigned char*) calloc(sizeof(unsigned char), filesize); fread(buffer, sizeof(unsigned char), filesize, mapfile); fseek(mapfile,0, SEEK_SET); MD5Update(&ctx, buffer, filesize); MD5Final(&ctx); memcpy(hash, ctx.digest, sizeof(ctx.digest)); free(buffer); buffer = NULL; } int get_function_idx(char* name){ int n, idx=-1; for (n=0; n> 6); } hash += (hash << 3); hash ^= (hash >> 11); hash += (hash << 15); return hash; } /* * Local variables: * c-indentation-style: bsd * c-basic-offset: 2 * indent-tabs-mode: nil * End: * * vim: tabstop=2 expandtab shiftwidth=2: */