Initial commit for the public release of the KYNrefrev XSPEC model.

lmodel.dat

+kynrefrev   36  0.         1.0e20         c_KYNrefrev   add  0
+a/M      GM/c     1.      0.      0.      1.      1.      0.2
+theta_o  deg      30.     0.      0.      89.     89.     5.
+rin      GM/c^2   1.      1.      1.      1000.   1000.   -0.5
+$ms      1.
+rout     GM/c^2   1000.   1.      1.      1000.   1000    -1.0
+phi      deg      0.      -180.   -180.   180.    180.    -1.
+dphi     deg      360.    0.      0.      360.    360.    -1.
+M/M8     " "      0.1     1e-8    1e-8    1e+3    1e+3    5.
+height   GM/c^2   3.      -20.    -20.    100.    100.    -1.
+PhoIndex " "      2.0     1.0     1.0     3.      3.      -0.1
+L/Ledd   " "      0.001   -1e+10  -1e+10  -1e-10  -1e-10  10.
+Np:Nr    " "      1.0     0.      0.      10.     10.     -0.1
+density  " "      1.      1e-8    1e-8    1e+8    1e+8    5.
+den_prof " "      0.      -5.     -5.     0.      0.      0.1
+abun     " "      1.      0.1     0.1     20.     20.     0.01
+alpha    GM/c^2   -6.     -100.   -100.   100.    100.    1.
+beta     GM/c^2   0.      -100.   -100.   100.    100.    1.
+rcloud   GM/c^2   0.      0.      0.      100.    100.    1.
+zshift   " "      0.      -0.999  -0.999  10.     10.     -0.1
+limb     " "      0.      0.      0.      2.      2.      -1.
+$tab     2.
+$sw      2.
+ntable   " "      80.     0.      0.      99.     99.     -1.
+nrad     " "      -4488.  -10000. -10000. 10000.  10000.  -100.
+division " "      -1.     -1.     -1.      1.      1.      -1.
+nphi     " "      180.    1.      1.      20000.  20000.  -100.
+deltaT   GM/c^3   1.      1e-3    1e-3    10.     10.     -0.05
+nt_ratio " "      1.      1.      1.      10.     10.     -1.
+t1/f1/E1 s/Hz/keV 0.3     -1e8    -1e8    1e16    1e16    -1.
+t2/f2/E2 s/Hz/keV 0.8     -1e8    -1e8    1e16    1e16    -1.
+Eref1    keV      1.      -1.     -1.     100.    100.    -1.
+Eref2    keV      3.      0.      0.      100.    100.    -1.
+dt/Af    s/-      0.      -1e8    -1e8    1e8     1e8     -1.
+k/qf     " "      1.      -1e8    -1e8    1e8     1e8     -1.
+$xsw     16.
+nthreads " "      4.      1.      1.      100.    100.    -1.

other/Makefile

+SRC_KYNREFREV=xside_threads.c xside.c libxspec.c fft_reverberation.c xskynrefrev.c
+
+CC=gcc
+INCLUDE=-I/usr/local/share/heasoft/x86_64-unknown-linux-gnu-libc2.19/include
+LIBRARY_PATH=-L/usr/local/share/heasoft/x86_64-unknown-linux-gnu-libc2.19/lib
+LIBRARY= -lcfitsio_3.37 -lpthread
+CFLAGS=-fPIC -O3 -Wall -DOUTSIDE_XSPEC -lm
+
+kynrefrev: $(SRC_KYNREFREV)
+		$(CC) $(CFLAGS) $(INCLUDE) -o $@ $(SRC_KYNREFREV) $(LIBRARY_PATH) $(LIBRARY)

other/libxspec.c

+#include <stdio.h>
+
+/*
+   Some routines defined in xspec and used in KY local models xs*.c
+   they are needed when we use routines outside xspec...
+*/
+
+// This function should return string with path to KBHtablesNN.fits
+// We define the string to be "./" (the tables are in working dir)
+char* FGMODF(void){
+  return("./");
+}
+
+// This function should return string with path to KBHtablesNN.fits
+// defined in KYDIR
+// We define the string to be ""
+char* FGMSTR(char* dname){
+  return("");
+}
+
+// This function should define and set some variable inside xspec
+void FPMSTR(const char* value1, const char* value2){
+  return;
+}
+
+// Next subroutine should write a text to the terminal
+int xs_write(char* text,int  idest){
+  fprintf(stdout,"%s\n",text);
+  return(0);
+}

xside_threads.c

+//#include <stdio.h>
+#include <unistd.h>
+#include <math.h>
+#include <pthread.h>
+#include <string.h>
+#include <stdlib.h>
+
+#define STAT_WAITING       1
+#define STAT_WORKING       2
+#define STAT_FINISHED      3
+#define STAT_NO_MORE       4
+
+struct t_params
+{
+	int    id;
+	int    status;	
+	int    iloop;
+	int    ne_loc;
+	int    ne;
+	long    nt;
+	double *fce;
+	double *qce;
+	double *uce;
+	double *vce;
+	void   (*emissivity)();
+	void   (*element)();
+};
+
+static pthread_cond_t  thread_work = PTHREAD_COND_INITIALIZER;
+static pthread_cond_t  main_work = PTHREAD_COND_INITIALIZER;
+static pthread_mutex_t lock = PTHREAD_MUTEX_INITIALIZER;
+
+static pthread_t       *th;
+static struct t_params *params;
+
+
+void *thread(void *arg){
+
+ struct t_params *my_params;
+
+ my_params = (struct t_params*)arg;
+
+ pthread_mutex_lock(&lock);
+ while(my_params->status != STAT_NO_MORE){
+  if(my_params->status == STAT_WORKING){
+   pthread_mutex_unlock(&lock);
+   my_params->element(my_params->iloop,my_params->ne_loc,my_params->ne,
+            my_params->nt,my_params->fce,my_params->qce,my_params->uce,
+            my_params->vce,my_params->emissivity);
+   pthread_mutex_lock(&lock);
+   my_params->status = STAT_FINISHED;
+   pthread_cond_signal(&main_work);
+  }else pthread_cond_wait(&thread_work, &lock);
+ };
+ pthread_mutex_unlock(&lock);
+ pthread_exit(NULL);
+}
+
+
+void create_ide_threads(const int nthreads, const int nthreads_old, 
+                        const int ne, const int ne_old,
+                        const long nt, const long nt_old, 
+                        const int polar, const int polar_old){
+
+int   i;
+void  *return_pointer;
+
+// Let's free pointers if we change nthreads, ne, nt or polar
+ if(nthreads_old > 1 && (nthreads != nthreads_old || ne != ne_old || 
+                         nt != nt_old || (polar == 0 && polar_old == 1))){
+  pthread_mutex_lock(&lock); 
+  if(nthreads != nthreads_old || ne != ne_old || nt != nt_old) 
+   for(i = 0; i < nthreads_old; i++) free((double *) params[i].fce);
+  if(polar_old) for(i = 0; i < nthreads_old; i++){
+    free((double *) params[i].qce);
+    free((double *) params[i].uce);
+    free((double *) params[i].vce);
+  };
+  pthread_mutex_unlock(&lock); 
+  if(nthreads != nthreads_old){
+   pthread_mutex_lock(&lock); 
+   for(i = 0; i < nthreads_old; i++) params[i].status = STAT_NO_MORE;
+   pthread_cond_broadcast(&thread_work);
+   pthread_mutex_unlock(&lock); 
+   for(i = 0; i < nthreads_old; i++) pthread_join(*(th + i), &return_pointer);
+   free(th);
+   free(params);
+  };
+ };
+// Let's create pointers if we change nthreads, ne, nt or polar
+ if(nthreads > 1){
+  if(nthreads != nthreads_old){
+   th = (pthread_t *)malloc(nthreads * sizeof(pthread_t));
+   params = (struct t_params *)malloc(nthreads * sizeof(struct t_params));
+   for(i = 0; i < nthreads; i++){
+    params[i].id = i;
+    params[i].status = STAT_WAITING;
+    pthread_create(th + i, NULL, thread, (void *)(params + i));
+   }; 
+  };
+  if(nthreads != nthreads_old || ne != ne_old || nt != nt_old)
+   for(i = 0; i < nthreads; i++){
+    params[i].ne = ne;
+    params[i].nt = nt;
+    params[i].fce = (double *)malloc(ne*nt*sizeof(double));
+   };
+  if(polar && (nthreads != nthreads_old || ne != ne_old || nt != nt_old || 
+     polar_old == 0)){
+   for(i = 0; i < nthreads; i++){
+    params[i].qce = (double *)malloc(ne*nt*sizeof(double));
+    params[i].uce = (double *)malloc(ne*nt*sizeof(double));
+    params[i].vce = (double *)malloc(ne*nt*sizeof(double));
+   };
+  };
+ };
+}
+
+
+void ide_threads(const int nthreads, const int nloops, const int polar, 
+                 const int ne_loc, const int ne, const long nt, double *fc, 
+                 double *qc, double *uc, double *vc, void (*emissivity)(), 
+                 void (*element)()){
+                  
+ int       i, ie;
+ long      it;
+ int       iloop, finished;
+ 
+ finished = 0;
+ iloop = 1;
+ pthread_mutex_lock(&lock);
+ for(i = 0; i < nthreads; i++){
+  params[i].status = STAT_WAITING;
+  params[i].ne_loc = ne_loc;
+  for(ie = 0; ie < ne; ie++) for(it = 0; it < nt; it++){
+   *(params[i].fce+ie+ne*it)=0.;
+   if(polar){
+    *(params[i].qce+ie+ne*it)=0.;
+    *(params[i].uce+ie+ne*it)=0.;
+    *(params[i].vce+ie+ne*it)=0.;
+   };
+  };
+  params[i].emissivity=emissivity;
+  params[i].element=element;
+ };
+ while(finished < nloops){
+  for(i = 0; (i < nthreads) && (iloop <= nloops); i++) 
+   if(params[i].status != STAT_WORKING){
+    params[i].status = STAT_WORKING;
+    params[i].iloop = iloop;
+    iloop++;
+  };
+  pthread_cond_broadcast(&thread_work);
+  pthread_cond_wait(&main_work, &lock);
+  for(i = 0; i < nthreads; i++) if(params[i].status == STAT_FINISHED){
+   finished++;
+   params[i].status = STAT_WAITING;
+  };
+ };
+ for(i = 0; i < nthreads; i++){
+  for(ie = 0; ie < ne; ie++) for(it = 0; it < nt; it++){
+   *(fc+ie+ne*it) += *(params[i].fce+ie+ne*it);
+   if(polar){
+    *(qc+ie+ne*it) += *(params[i].qce+ie+ne*it);
+    *(uc+ie+ne*it) += *(params[i].uce+ie+ne*it);
+    *(vc+ie+ne*it) += *(params[i].vce+ie+ne*it);
+   };
+  };
+ };
+ pthread_mutex_unlock(&lock);
+ return;
+}