util.c

   int j;
   long k;
   static long idum2=123456789;
   static long iy=0;
   static long iv[NTAB];
   double temp;
   
   if (*idum <= 0) {                 // Initialize.
       if (-(*idum) < 1) *idum=1;     // Be sure to prevent idum = 0.
       else *idum = -(*idum);
       idum2=(*idum);
       for (j=NTAB+7;j>=0;j--) {      // Load the shu.e table (after 8 warm-ups).
           k=(*idum)/IQ1;
           *idum=IA1*(*idum-k*IQ1)-k*IR1;
           if (*idum < 0) *idum += IM1;
           if (j < NTAB) iv[j] = *idum;
       }
       iy=iv[0];
   }
   k=(*idum)/IQ1;                    // Start here when not initializing.
   *idum=IA1*(*idum-k*IQ1)-k*IR1;    // Compute idum=(IA1*idum) % IM1 without
                                     // over ows by Schrage's method.
   if (*idum < 0) *idum += IM1;
   k=idum2/IQ2;
   idum2=IA2*(idum2-k*IQ2)-k*IR2;    // Compute idum2=(IA2*idum) % IM2 likewise.
   if (idum2 < 0) idum2 += IM2;
   j=iy/NDIV;                        // Will be in the range 0..NTAB-1.
   iy=iv[j]-idum2;                   // Here idum is shu.ed, idum and idum2 are
                                     // combined to generate output.
   iv[j] = *idum;
   if (iy < 1) iy += IMM1;
   if ((temp=AM*iy) > RNMX) return RNMX;  // Because users don't expect endpoint values.
   else return temp;
}


void ByteSwap(unsigned char * b, int n)
{
   register int i = 0;
   register int j = n-1;
   char temp;
   while (i<j)
   {
 		 temp = b[i];
         b[i] = b[j];
 		 b[j] = temp;
                 /* std::swap(b[i], b[j]); */
      i++, j--;
   }
}

int invMap(int nAstroPSolved,int *inv){
//inv[0] Row of Long
//inv[1] Row of Lat
//inv[2] Row of MuLong
//inv[3] Row of MuLat
    
    for(int i=0;i<4;i++)
        inv[i]=-1;
    if(nAstroPSolved<1 || nAstroPSolved>5)
        return 1;
    if(nAstroPSolved==2){
        inv[0]=0;
        inv[1]=1;
    }
    if(nAstroPSolved==3){
        inv[0]=1;
        inv[1]=2;
    }
    if(nAstroPSolved==4){
        inv[0]=0;
        inv[1]=1;
        inv[2]=2;
        inv[3]=3;
    }
    if(nAstroPSolved==5){
        inv[0]=1;
        inv[1]=2;
        inv[2]=3;
        inv[3]=4;
    }
    return 0;
}

void writeBinFiles(double* systemMatrix,long* matrixIndex,int* instrCol,double* knownTerms,char* wrfileDir,char* wpath, struct comData comlsqr, int debugMode){
    
    int nproc,myid;
    FILE *fpSM,*fpMI,*fpII,*fpKT,*fpBar;
    int nparam;
    
    MPI_Status status;
    char actpath[1024];
    size_t sizePath=1020;
    int extConstraint;
    int barConstraint;

    
    MPI_Comm_size(MPI_COMM_WORLD, &nproc);
    MPI_Comm_rank(MPI_COMM_WORLD, &myid);
    
    nparam=comlsqr.nAstroPSolved+comlsqr.nAttP+comlsqr.nInstrPSolved+comlsqr.nGlobP;
    getcwd(actpath,sizePath);
    chdir(wpath);
    if(!(chdir(wrfileDir)==0))  {
        printf("wrfile directory does not exist. Aborting\n");
        MPI_Abort(MPI_COMM_WORLD,1);
        exit(EXIT_FAILURE);
    }

    
    long int rowCounter=1;
    long int rowWritten=0;
    long int stillOpen=-1;
    int nOfFile=1;  // file numerati da 1 in avanti. Il corrente file da scrivere è nOfFile
    int lastStargsrId=-1;
    rowWritten=0;
	rowCounter=0;
	long rowCounterPrev=0;
    extConstraint=comlsqr.extConstraint;
    barConstraint=comlsqr.barConstraint;
	
    if(myid>0){
        MPI_Recv(&lastStargsrId, 1, MPI_INT, myid-1, 0, MPI_COMM_WORLD, &status);
        MPI_Recv(&rowCounterPrev, 1, MPI_LONG, myid-1, 1, MPI_COMM_WORLD, &status); //righe scritte nell'ultimo file
        MPI_Recv(&nOfFile, 1, MPI_INT, myid-1, 2, MPI_COMM_WORLD, &status);
//        printf("TP0 PE=%d ricevo da %d  lastStargsrId=%d rowCounterPrev=%ld nOfFile=%ld\n",myid,myid-1,lastStargsrId,rowCounterPrev,nOfFile);
    }
    if(lastStargsrId==1000*nOfFile-1)  // sono in fase di chiusura file
    {
        if(lastStargsrId!=matrixIndex[0]/comlsqr.nAstroPSolved ){ //rinomino il file precedente era alla fine, io sono in stella successiva!
            char filenameCoeffBinMatrixIndex[512];
            strcpy(filenameCoeffBinMatrixIndex,"Gsr_dpccu3dpctavugsrgsrsystemrow_0000_");
            char tmp1[20],tmp2[20],tmp3[20];
            sprintf(tmp1,"%09d",1000*(nOfFile-1));
            sprintf(tmp2,"%09d",lastStargsrId);
            sprintf(tmp3,"%07ld",rowCounterPrev);
            rowCounterPrev=0;
            
            strcat(filenameCoeffBinMatrixIndex, tmp1);
            strcat(filenameCoeffBinMatrixIndex, "_");
            strcat(filenameCoeffBinMatrixIndex, tmp2);
            strcat(filenameCoeffBinMatrixIndex, "_000000_nrows-");
            strcat(filenameCoeffBinMatrixIndex, tmp3);
            strcat(filenameCoeffBinMatrixIndex, "_MI.bin");

            char filenameCoeffBinSystemMatrix[512];
            strcpy(filenameCoeffBinSystemMatrix,"Gsr_dpccu3dpctavugsrgsrsystemrow_0000_");
            strcat(filenameCoeffBinSystemMatrix, tmp1);
            strcat(filenameCoeffBinSystemMatrix, "_");
            strcat(filenameCoeffBinSystemMatrix, tmp2);
            strcat(filenameCoeffBinSystemMatrix, "_000000_nrows-");
            strcat(filenameCoeffBinSystemMatrix, tmp3);
            strcat(filenameCoeffBinSystemMatrix, "_SM.bin");
            
            char filenameCoeffBinIstrIndex[512];
            strcpy(filenameCoeffBinIstrIndex,"Gsr_dpccu3dpctavugsrgsrsystemrow_0000_");
            strcat(filenameCoeffBinIstrIndex, tmp1);
            strcat(filenameCoeffBinIstrIndex, "_");
            strcat(filenameCoeffBinIstrIndex, tmp2);
            strcat(filenameCoeffBinIstrIndex, "_000000_nrows-");
            strcat(filenameCoeffBinIstrIndex, tmp3);
            strcat(filenameCoeffBinIstrIndex, "_II.bin");
            
            char filenameCoeffBinKnownTerms[512];
            strcpy(filenameCoeffBinKnownTerms,"Gsr_dpccu3dpctavugsrgsrsystemrow_0000_");
            strcat(filenameCoeffBinKnownTerms, tmp1);
            strcat(filenameCoeffBinKnownTerms, "_");
            strcat(filenameCoeffBinKnownTerms, tmp2);
            strcat(filenameCoeffBinKnownTerms, "_000000_nrows-");
            strcat(filenameCoeffBinKnownTerms, tmp3);
            strcat(filenameCoeffBinKnownTerms, "_KT.bin");
            stillOpen=0;
            rename("MI.bin", filenameCoeffBinMatrixIndex);
            rename("SM.bin", filenameCoeffBinSystemMatrix);
            rename("II.bin", filenameCoeffBinIstrIndex);
            rename("KT.bin", filenameCoeffBinKnownTerms);
            nOfFile++;
//            printf("TP3 PE=%d rinominati file %s %s %s %s\n",myid,filenameCoeffBinMatrixIndex, filenameCoeffBinSystemMatrix,filenameCoeffBinIstrIndex,filenameCoeffBinKnownTerms);

        }  // if(lastStargsrId!=matrixIndex[0]/comlsqr.nAstroPSolved
    } // if(lastStargsrId==1000*nOfFile-1
 
       
    fpMI=fopen("MI.bin","ab");
    fpSM=fopen("SM.bin","ab");
    fpII=fopen("II.bin","ab");
    fpKT=fopen("KT.bin","ab");
//    printf("TP4 PE=%d aperti in append  files MI.bin ecc\n",myid);
    for(long int i=0;i<comlsqr.mapNoss[myid];i++){
       if(i==comlsqr.mapNoss[myid]-1 ){ //  sono all'ultima osservazione dell'attuale pe
//           printf("TP5 PE=%d  i=%ld\n",myid,i);
			rowCounter++;
			fwrite(&matrixIndex[rowWritten*comlsqr.multMI],sizeof(long int),rowCounter*comlsqr.multMI,fpMI);
                     
			for(int q=0;q<nparam;q++) {
					fwrite(&systemMatrix[rowWritten*nparam+rowCounter*q],sizeof(double),rowCounter,fpSM);
			}
           
			if(comlsqr.nInstrPSolved>0) fwrite(&instrCol[rowWritten*comlsqr.nInstrPSolved],sizeof(int),rowCounter*comlsqr.nInstrPSolved,fpII);
					
			fwrite(&knownTerms[rowWritten],sizeof(double),rowCounter,fpKT);
			fclose(fpMI);
			fclose(fpSM);
			fclose(fpII);
			fclose(fpKT);
//            printf("TP6 PE=%d  i=%ld, rowCounter=%ld rowWritten=%ld rowCounterPrev=%ld\n",myid,i,rowCounter,rowWritten,rowCounterPrev);

            if(myid==nproc-1){
                 lastStargsrId=matrixIndex[comlsqr.mapNoss[myid]*2-comlsqr.multMI]/comlsqr.nAstroPSolved;
//                printf("TP7 PE=%d  i=%ld lastStargsrId=%d\n",myid,i,lastStargsrId);
               char filenameCoeffBinMatrixIndex[512];
               strcpy(filenameCoeffBinMatrixIndex,"Gsr_dpccu3dpctavugsrgsrsystemrow_0000_");
               char tmp1[20],tmp2[20],tmp3[20];
               sprintf(tmp1,"%09d",1000*(nOfFile-1));
               sprintf(tmp2,"%09d",lastStargsrId);
               sprintf(tmp3,"%07ld",rowCounter+rowCounterPrev);
               strcat(filenameCoeffBinMatrixIndex, tmp1);
               strcat(filenameCoeffBinMatrixIndex, "_");
               strcat(filenameCoeffBinMatrixIndex, tmp2);
               strcat(filenameCoeffBinMatrixIndex, "_000000_nrows-");
               strcat(filenameCoeffBinMatrixIndex, tmp3);
               strcat(filenameCoeffBinMatrixIndex, "_MI.bin");

               char filenameCoeffBinSystemMatrix[512];
               strcpy(filenameCoeffBinSystemMatrix,"Gsr_dpccu3dpctavugsrgsrsystemrow_0000_");
               strcat(filenameCoeffBinSystemMatrix, tmp1);
               strcat(filenameCoeffBinSystemMatrix, "_");
               strcat(filenameCoeffBinSystemMatrix, tmp2);
               strcat(filenameCoeffBinSystemMatrix, "_000000_nrows-");
               strcat(filenameCoeffBinSystemMatrix, tmp3);
               strcat(filenameCoeffBinSystemMatrix, "_SM.bin");
               char filenameCoeffBinIstrIndex[512];
               strcpy(filenameCoeffBinIstrIndex,"Gsr_dpccu3dpctavugsrgsrsystemrow_0000_");
               strcat(filenameCoeffBinIstrIndex, tmp1);
               strcat(filenameCoeffBinIstrIndex, "_");
               strcat(filenameCoeffBinIstrIndex, tmp2);
               strcat(filenameCoeffBinIstrIndex, "_000000_nrows-");
               strcat(filenameCoeffBinIstrIndex, tmp3);
               strcat(filenameCoeffBinIstrIndex, "_II.bin");
               char filenameCoeffBinKnownTerms[512];
               strcpy(filenameCoeffBinKnownTerms,"Gsr_dpccu3dpctavugsrgsrsystemrow_0000_");
               strcat(filenameCoeffBinKnownTerms, tmp1);
               strcat(filenameCoeffBinKnownTerms, "_");
               strcat(filenameCoeffBinKnownTerms, tmp2);
               strcat(filenameCoeffBinKnownTerms, "_000000_nrows-");
               strcat(filenameCoeffBinKnownTerms, tmp3);
               strcat(filenameCoeffBinKnownTerms, "_KT.bin");
               rename("MI.bin", filenameCoeffBinMatrixIndex);
               rename("SM.bin", filenameCoeffBinSystemMatrix);
               rename("II.bin", filenameCoeffBinIstrIndex);
               rename("KT.bin", filenameCoeffBinKnownTerms);
                rowCounterPrev=rowCounter;  //inutile
//               printf("TP7.1 PE=%d  i=%ld, rowCounter=%ld rowWritten=%ld  rinominati i files MI.bin ecc FINITO!!\n",myid,i,rowCounter,rowWritten);
                break; 
           } //if(myid==nproc-1)
           rowCounterPrev+=rowCounter;
//           printf("TP7.2 PE=%d  i=%ld, rowCounter=%ld rowWritten=%ld rowCounterPrev=%ld PROX PE!!\n",myid,i,rowCounter,rowWritten,rowCounterPrev);
           break;  //esce dal ciclo for
			
        }//if(i==comlsqr.mapNoss[myid]-1 )
			
       if ( (matrixIndex[i*comlsqr.multMI]/comlsqr.nAstroPSolved) < 1000*nOfFile ) {
					rowCounter++;
					continue;
		} else
        {
            fwrite(&matrixIndex[rowWritten*comlsqr.multMI],sizeof(long int),rowCounter*comlsqr.multMI,fpMI);
            fclose(fpMI);
            
            for(int q=0;q<nparam;q++){
                fwrite(&systemMatrix[rowWritten*nparam+rowCounter*q],sizeof(double),rowCounter,fpSM);
            }

            fclose(fpSM);
            
            if(comlsqr.nInstrPSolved>0) fwrite(&instrCol[rowWritten*comlsqr.nInstrPSolved],sizeof(int),rowCounter*comlsqr.nInstrPSolved,fpII);
            fclose(fpII);
            
            fwrite(&knownTerms[rowWritten],sizeof(double),rowCounter,fpKT);
            fclose(fpKT);
            char filenameCoeffBinMatrixIndex[512];
			char filenameCoeffBinSystemMatrix[512];
			char filenameCoeffBinInstrIndex[512];
			char filenameCoeffBinKnownTerms[512];
			strcpy(filenameCoeffBinMatrixIndex,"Gsr_dpccu3dpctavugsrgsrsystemrow_0000_");
			strcpy(filenameCoeffBinSystemMatrix,"Gsr_dpccu3dpctavugsrgsrsystemrow_0000_");
			strcpy(filenameCoeffBinInstrIndex,"Gsr_dpccu3dpctavugsrgsrsystemrow_0000_");
			strcpy(filenameCoeffBinKnownTerms,"Gsr_dpccu3dpctavugsrgsrsystemrow_0000_");
			
            char tmp1[20],tmp2[20],tmp3[20];
			sprintf(tmp1,"%09d",1000*(nOfFile-1));
			
            sprintf(tmp2,"%09d",1000*nOfFile-1);
			sprintf(tmp3,"%07ld",rowCounter+rowCounterPrev);
			strcat(filenameCoeffBinMatrixIndex, tmp1);
			strcat(filenameCoeffBinMatrixIndex, "_");
			strcat(filenameCoeffBinMatrixIndex, tmp2);
			strcat(filenameCoeffBinMatrixIndex, "_000000_nrows-");
			strcat(filenameCoeffBinMatrixIndex, tmp3);
			strcat(filenameCoeffBinMatrixIndex, "_MI.bin");

			strcat(filenameCoeffBinSystemMatrix, tmp1);
			strcat(filenameCoeffBinSystemMatrix, "_");
			strcat(filenameCoeffBinSystemMatrix, tmp2);
			strcat(filenameCoeffBinSystemMatrix, "_000000_nrows-");
			strcat(filenameCoeffBinSystemMatrix, tmp3);
			strcat(filenameCoeffBinSystemMatrix, "_SM.bin");
			strcat(filenameCoeffBinInstrIndex, tmp1);
			strcat(filenameCoeffBinInstrIndex, "_");
			strcat(filenameCoeffBinInstrIndex, tmp2);
			strcat(filenameCoeffBinInstrIndex, "_000000_nrows-");
			strcat(filenameCoeffBinInstrIndex, tmp3);
			strcat(filenameCoeffBinInstrIndex, "_II.bin");
			strcat(filenameCoeffBinKnownTerms, tmp1);
            strcat(filenameCoeffBinKnownTerms, "_");
            strcat(filenameCoeffBinKnownTerms, tmp2);
            strcat(filenameCoeffBinKnownTerms, "_000000_nrows-");
            strcat(filenameCoeffBinKnownTerms, tmp3);
            strcat(filenameCoeffBinKnownTerms, "_KT.bin");
            rename("MI.bin", filenameCoeffBinMatrixIndex);
            rename("SM.bin", filenameCoeffBinSystemMatrix);
            rename("II.bin", filenameCoeffBinInstrIndex);
            rename("KT.bin", filenameCoeffBinKnownTerms);
//            printf("TP9 PE=%d  i=%ld Rinominati sui file GSR_... MI.bin ecc %s \n",myid,i,filenameCoeffBinMatrixIndex );
            fpMI=fopen("MI.bin","wb");
            fpSM=fopen("SM.bin","wb");
            fpII=fopen("II.bin","wb");
            fpKT=fopen("KT.bin","wb");
            rowWritten+=rowCounter;
			rowCounter=1;  // sono già sulla successiva stella
//            printf("TP10 PE=%d aperti in wb  files MI.bin ecc\n",myid);
            rowCounterPrev=0;
			
            nOfFile++;
        }//else //Scrivere su file
    }
    if(myid<nproc-1){
        lastStargsrId=matrixIndex[comlsqr.mapNoss[myid]*2-comlsqr.multMI]/comlsqr.nAstroPSolved;
        
        MPI_Send(&lastStargsrId, 1, MPI_INT, myid+1, 0, MPI_COMM_WORLD);
        MPI_Send(&rowCounterPrev, 1, MPI_LONG, myid+1, 1, MPI_COMM_WORLD);
        MPI_Send(&nOfFile, 1, MPI_INT, myid+1, 2, MPI_COMM_WORLD);


    }

    if(extConstraint && myid==0){
        /// write extConstr for stars
        long nStar=comlsqr.nStar;
        int nEqExtConstr=comlsqr.nEqExtConstr;
        long nDegFreedomAtt=comlsqr.nDegFreedomAtt;
        int startStar=0;
        int endStar=0;
        int numOfExtStarinFile=1000000;
        int numOfExtStartoWrite;
        double randVal;
        short nAstroPSolved=comlsqr.nAstroPSolved;
        for(int nwr=0;nwr<nStar;nwr+=numOfExtStarinFile){   
           int localCounter=0;
           if(nStar>numOfExtStarinFile+nwr)numOfExtStartoWrite=numOfExtStarinFile;
            else numOfExtStartoWrite=nStar-nwr;
            startStar=nwr;
            endStar=startStar+numOfExtStartoWrite-1;
            double *buffArray;
            buffArray=(double *)calloc(numOfExtStartoWrite*nAstroPSolved*nEqExtConstr,sizeof(double));
 
            for(int j=0;j<nEqExtConstr;j++)
                for(int i=0;i<numOfExtStartoWrite*nAstroPSolved;i++){
                    randVal=(((double)rand())/RAND_MAX)*2 - 1.0;
                    if(nAstroPSolved==3 && i%nAstroPSolved==0) randVal=0.;
                    if(nAstroPSolved==4 && i%nAstroPSolved>=2) randVal=0.;
                    if(nAstroPSolved==5 && i%nAstroPSolved==0) randVal=0.;
                    if(nAstroPSolved==5 && i%nAstroPSolved>2 && j<3) randVal=0.;
                    buffArray[localCounter]=randVal;
                    localCounter++;
                }
            
            
            char filenameExtConstrAstro[512];
            char tmp1[20],tmp2[20],tmp3[20];
            sprintf(tmp1,"%09d",startStar);
            sprintf(tmp2,"%09d",endStar);
            sprintf(tmp3,"%09d",nEqExtConstr);
            strcpy(filenameExtConstrAstro,"Gsr_nullspaceastrofit_0000_");
            strcat(filenameExtConstrAstro, tmp1);
            strcat(filenameExtConstrAstro, "_");
            strcat(filenameExtConstrAstro, tmp2);
            strcat(filenameExtConstrAstro, "_000000_nrows-");
            strcat(filenameExtConstrAstro, tmp3);
            strcat(filenameExtConstrAstro,".bin");
            fpSM=fopen(filenameExtConstrAstro,"wb");
            fwrite(buffArray,sizeof(double),nEqExtConstr*numOfExtStartoWrite*nAstroPSolved,fpSM);
            fclose(fpSM);
            free(buffArray);
        } // for(in nwr=..
        /// write extConstr for Att
        double *buffArray;
        buffArray=(double *)calloc(nDegFreedomAtt,sizeof(double));
        for(int i=0;i<nDegFreedomAtt;i++){
            randVal=(((double)rand())/RAND_MAX)*2 - 1.0;
            buffArray[i]=randVal;
        }
        char filenameExtConstrAtt[512];
        strcpy(filenameExtConstrAtt,"Gsr_nullspaceattitudefit_0000_000000_nrows-0660601_EC.bin");
        fpSM=fopen(filenameExtConstrAtt,"wb");
        fwrite(buffArray,sizeof(double),nDegFreedomAtt,fpSM);
        fclose(fpSM);
        free(buffArray);

    }//if (extConstraint

    if(barConstraint && myid==0){
        /// write extConstr for stars
        long nStar=comlsqr.nStar;
        int nEqBarConstr=comlsqr.nEqBarConstr;
        int startStar=0;
        int endStar=0;
        int numOfBarStarinFile=1000000;
        int numOfBarStartoWrite;
        double randVal;
        short nAstroPSolved=comlsqr.nAstroPSolved;
        for(int nwr=0;nwr<nStar;nwr+=numOfBarStarinFile){
            int localCounter=0;
            if(nStar>numOfBarStarinFile+nwr)numOfBarStartoWrite=numOfBarStarinFile;
            else numOfBarStartoWrite=nStar-nwr;
            startStar=nwr;
            endStar=startStar+numOfBarStartoWrite-1;
            double *buffArray;
            buffArray=(double *)calloc(numOfBarStartoWrite*nAstroPSolved*nEqBarConstr,sizeof(double));
            
            for(int j=0;j<nEqBarConstr;j++)
            for(int i=0;i<numOfBarStartoWrite*nAstroPSolved;i++){
                randVal=(((double)rand())/RAND_MAX)*2 - 1.0;
                if(nAstroPSolved==3 && i%nAstroPSolved==0) randVal=0.;
                if(nAstroPSolved==4 && i%nAstroPSolved>=2) randVal=0.;
                if(nAstroPSolved==5 && i%nAstroPSolved==0) randVal=0.;
                if(nAstroPSolved==5 && i%nAstroPSolved>2 && j<3) randVal=0.;
                buffArray[localCounter]=randVal;
                localCounter++;
            }
            
            
            char filenameBarConstrAstro[512];
            char tmp1[20],tmp2[20],tmp3[20];
            sprintf(tmp1,"%09d",startStar);
            sprintf(tmp2,"%09d",endStar);
            sprintf(tmp3,"%09d",nEqBarConstr);
            strcpy(filenameBarConstrAstro,"Gsr_barconstrastrofit_0000_");
            strcat(filenameBarConstrAstro, tmp1);
            strcat(filenameBarConstrAstro, "_");
            strcat(filenameBarConstrAstro, tmp2);
            strcat(filenameBarConstrAstro, "_000000_nrows-");
            strcat(filenameBarConstrAstro, tmp3);
            strcat(filenameBarConstrAstro,".bin");
            fpBar=fopen(filenameBarConstrAstro,"wb");
            fwrite(buffArray,sizeof(double),nEqBarConstr*numOfBarStartoWrite*nAstroPSolved,fpBar);
            fclose(fpBar);
            free(buffArray);
        } // for(in nwr=..
        
    }//if (barConstraint

    MPI_Barrier(MPI_COMM_WORLD);
    chdir(wpath);
    chdir(actpath);
}

int cmpfunc (const void * a, const void * b)
{
    return ( *(int*)a - *(int*)b );
}

int randint(int max) {
    return (int) (max*(rand()/(RAND_MAX+1.0)));
}

// This function generates a psuedo-random integer n, where min <= n <= max
int randint1(int min, int max) {
    if (min>max) {
        return max+(int)((min-max+1)*(rand()/(RAND_MAX+1.0)));
    } else {
        return min+(int)((max-min+1)*(rand()/(RAND_MAX+1.0)));
    }
}

// This function generates a psuedo-random long integer n, where 0 <= n < max
long randlong(long max) {
    return (long) (max*(rand()/(RAND_MAX+1.0)));
}

// This function generates a psuedo-random long integer n, where min <= n <= max
long randlong1(long min, long max) {
    if (min>max) {
        return max+(long)((min-max+1)*(rand()/(RAND_MAX+1.0)));
    } else {
        return min+(long)((max-min+1)*(rand()/(RAND_MAX+1.0)));
    }
}

// This function computes the instrument hash given the number of FoV, CCD, PixelColumn and TimeInterval.
// Since these numbers must fill in 1, 8, 11, and 11 bits respectively, they are checked against appropriate
// maximum values, and if their values are greater, the function returns -1.
// If the function is called with FoV, CCD, PixelColumn and TimeInterval it returns the instrId hash,
// if it is called with nFoVs, nCCDs, nPixelColumns and nTimeIntervals it returns the instrSetUp hash.
// NB: the FoV is treated differently since it can be 1 or 2, so what is actually hashed is nFoVs-1 to keep it into 1 bit.
long instr_hash(int FoV, int CCD, int PixelColumn, int TimeInterval) {
    short CCDOffset = 1;
    short PixelOffset = 9;
    short TimeOffset = 20;
    short FoVMaxValue = ((short) pow(2.0, CCDOffset));
    short CCDMaxValue = ((short) pow(2.0, PixelOffset-CCDOffset));
    short PixelMaxValue = ((short) pow(2.0, TimeOffset-PixelOffset));
    short TimeMaxValue = ((short) pow(2.0, TimeOffset-PixelOffset));
    
    if(FoV>FoVMaxValue || CCD>CCDMaxValue || PixelColumn>PixelMaxValue || TimeInterval>TimeMaxValue)
        return -1; // Some parameter is greater than its maximum possible value. The hash cannot be calculated.
    else
        return (((long) FoV-1)
                | (((long) CCD) << (CCDOffset))
                | (((long) PixelColumn) << (PixelOffset))
                | (((long) TimeInterval) << (TimeOffset)));
}


// This function extracts an integer in the range pos_min..pos_max which points to a position of the array <values>.
// This implies that values is an array with pos_max+1 positions.
void swap(long *i, long *j) {
    long t = *i;
    *i = *j;
    *j = t;
}

int fill_extract(long *values, long *pos_min, long pos_max, long *number) {
    long position;
    
    if(*pos_min==pos_max-1)
        position=*pos_min;
    else
        position=randlong1(*pos_min, pos_max-1);
    *number=values[position];
    swap(&values[*pos_min], &values[position]);
    (*pos_min)++;
    if(*pos_min>=pos_max) 
        return 1;
    return 0;
}

struct nullSpace cknullSpace(double * systemMatrix,long * matrixIndex,double *attNS,struct comData  comlsqr){
    struct nullSpace results;
    int nproc, myid;
    long nunkSplitNS;
    double * nullSpaceVect;
    double * prodNS;
    long nElements, nStar;
    int nEqExtConstr;
    int nparam;
    int nLocalStar;
    int firstStarConstr,lastStarConstr;
    int nOfElextObs,numOfExtStar,numOfExtAttCol,startingAttColExtConstr;
    short nAstroPSolved,nAttParAxis,nAttAxes;
    double *nullSpaceFPN;
    double *productNS;
    int npeSend,npeRecv;
    double sum,extConstrW;
    long sumVer;
    int setBound[4];
    long int l1,k,l,j;
    long int nDegFreedomAtt,localAstroMax,offsetAttParam;
    long int *mapNoss, *mapNcoeff;
    time_t seconds[2], tot_sec;

    int **mapStar;
    
    MPI_Status status;
    MPI_Request req1;

    MPI_Comm_size(MPI_COMM_WORLD, &nproc);
    MPI_Comm_rank(MPI_COMM_WORLD, &myid);
    
    
    nEqExtConstr=comlsqr.nEqExtConstr;
    firstStarConstr=comlsqr.firstStarConstr;
    lastStarConstr=comlsqr.lastStarConstr;
    nAstroPSolved=comlsqr.nAstroPSolved;
    nOfElextObs=comlsqr.nOfElextObs;
    nDegFreedomAtt=comlsqr.nDegFreedomAtt;
    nAttParAxis=comlsqr.nAttParAxis;
    localAstroMax=comlsqr.VrIdAstroPDimMax*nAstroPSolved;
    offsetAttParam=comlsqr.offsetAttParam;
    numOfExtStar=comlsqr.numOfExtStar;
    nAttAxes=comlsqr.nAttAxes;
    mapNcoeff=comlsqr.mapNcoeff;
    mapNoss=comlsqr.mapNoss;
    numOfExtAttCol=comlsqr.numOfExtAttCol;
    startingAttColExtConstr=comlsqr.startingAttColExtConstr;
    setBound[0]=comlsqr.setBound[0];
    setBound[1]=comlsqr.setBound[1];
    setBound[2]=comlsqr.setBound[2];
    setBound[3]=comlsqr.setBound[3];
    extConstrW=comlsqr.extConstrW;
    int multMI=comlsqr.multMI;
    mapStar=comlsqr.mapStar;
    nStar=comlsqr.nStar;
    int nAttParam=comlsqr.nAttParam;
    int nInstrParam=comlsqr.nInstrParam;
    int nGlobalParam=comlsqr.nGlobalParam;
    
// errore ==>    nparam=nAstroPSolved+comlsqr.nAttP+comlsqr.nInstrP+comlsqr.nGlobP;
    nparam=nAstroPSolved+comlsqr.nAttP+comlsqr.nInstrPSolved+comlsqr.nGlobP;
    nLocalStar=comlsqr.mapStar[myid][1]-comlsqr.mapStar[myid][0]+1;
    
    nElements = mapNoss[myid]+nEqExtConstr;
    
    nullSpaceFPN = (double *) calloc(nAstroPSolved*nEqExtConstr, sizeof(double));
    if (!nullSpaceFPN)
        exit(err_malloc("nullSpaceFPN",myid));
     
    productNS = (double *) calloc(nElements, sizeof(double));
    if (!productNS)
        exit(err_malloc("productNS",myid));

     
       
    npeSend=myid-1;
    npeRecv=myid+1;
        
    for(int i=0; i<nEqExtConstr;i++){
        if(myid>0){
        MPI_Isend(&systemMatrix[mapNoss[myid]*nparam+i*nOfElextObs],nAstroPSolved, MPI_DOUBLE, npeSend, 1,MPI_COMM_WORLD, &req1);
        }
        if(myid<nproc-1){
            MPI_Recv(&nullSpaceFPN[nAstroPSolved*i], nAstroPSolved, MPI_DOUBLE, npeRecv, 1,MPI_COMM_WORLD, &status);
        }
         if(myid>0) MPI_Wait(&req1,&status);
        
        MPI_Barrier(MPI_COMM_WORLD);
    }
        
        
        
    prodNS = (double *) calloc(nElements, sizeof(double));
    if (!prodNS)
            exit(err_malloc("prodNS",myid));

    
    
    nunkSplitNS=localAstroMax + nAttParam+nInstrParam+nGlobalParam;
    
   
    nullSpaceVect = (double *) calloc(nunkSplitNS, sizeof(double));
    if (!nullSpaceVect) exit(err_malloc("nullSpaceVect",myid));
    
    for(int ic=0; ic<nEqExtConstr;ic++){
        seconds[0]=time(NULL);
        for(int j1=localAstroMax;j1<nunkSplitNS;j1++)
            nullSpaceVect[j1]=0.0;
            for(int j1=0;j1<nElements;j1++)
                productNS[j1]=0.0;
        for(int j1=0;j1<(lastStarConstr-firstStarConstr+1)*nAstroPSolved;j1++){
            nullSpaceVect[j1]= systemMatrix[mapNoss[myid]*nparam+j1+ic*comlsqr.nOfElextObs]/extConstrW;
        }
        if(comlsqr.mapStar[myid][1]>lastStarConstr){
            for(int j=0;j<nAstroPSolved;j++)
                nullSpaceVect[(lastStarConstr-firstStarConstr+1)*nAstroPSolved+j]=nullSpaceFPN[ic*nAstroPSolved+j]/extConstrW;
        }

        if(ic<3)
        {
            for (int m=0;m<nDegFreedomAtt;m++)
                    nullSpaceVect[localAstroMax+ic*nDegFreedomAtt+m]=comlsqr.nullSpaceAttfact/extConstrW;
        } else{
            for (int m=0;m<nDegFreedomAtt;m++)
                nullSpaceVect[localAstroMax+(ic-3)*nDegFreedomAtt+m]=attNS[m]/extConstrW;
            
        }
////#pragma omp parallel private(myid, sum, sumVer  k, l1, l2, l, j,tid,nthreads) shared(mapNoss,comlsqr,nullSpaceVect,systemMatrix,matrixIndex,productNS)
{
////#ifdef OMP
///            tid = omp_get_thread_num();
///            nthreads = omp_get_num_threads();
//            if(myid==39) printf("TP00.1 PE=%d tid=%d nthreads=%d mapNoss[myid]=%ld, nElements=%ld, multMI=%d\n",myid,tid,nthreads,mapNoss[myid],nElements,multMI);
//            if(myid==39) printf("TP00.1 PE=%d  mapNoss[myid]=%ld\n",myid,mapNoss[myid]);
////#endif
///#pragma omp for
        for(int i=0;i<mapNoss[myid];i++){
                sum=0.0;
                sumVer=0;
                k=0;
                l1=nparam*i;
                int lset=0;

            // sumAstro=(nAstroPSolved*nAstroPSolved*((j-1)*nStar+(multMI*i)/nAstroPSolved)+nAstroPSolved*(nAstroPSolved+1)/2)
            // sumAtt  =4*((j-1)*nDegFreedomAtt+(multMI*i+1))+6
            // sumVer  = sumAstro + sumAtt

            long chkSumVer=(nAstroPSolved*nAstroPSolved*(ic*nStar+matrixIndex[multMI*i]/nAstroPSolved)+(nAstroPSolved*(nAstroPSolved-1))/2)+4*(ic*nDegFreedomAtt+matrixIndex[multMI*i+1]-nStar*nAstroPSolved)+6;
            
            
                for(l=l1;l<l1+setBound[2];l++){
                        if(lset<setBound[1])
                        {
                            if(lset==0)
                            {
/*                                 if(multMI*i>=mapNoss[myid]*multMI){
                                    printf("ERROR: PE=%d i=%d lset=%d multMI*i=%d greather than mapNoss[myid]*multMI=%ld\n",myid,i,lset,multMI*i,mapNoss[myid]*multMI);
                                    continue;
			        }*/
                               long numOfStarPos=matrixIndex[multMI*i]/nAstroPSolved;
                                j=(numOfStarPos-comlsqr.mapStar[myid][0])*nAstroPSolved;
                            }
                            else
                                j++;
                        }
                        if(lset>=setBound[1] && lset<setBound[2])
                        {
                            if(((lset-setBound[1]) % nAttParAxis)==0) {
/*                                if(multMI*i>=mapNoss[myid]*multMI){
                                    printf("ERROR: PE=%d i=%d lset=%d multMI*i+(multMI-1)=%d greather than mapNoss[myid]*multMI=%ld\n",myid,i,lset,multMI*i+(multMI-1),mapNoss[myid]*multMI);
                                    continue;
                                }*/
                                j=matrixIndex[multMI*i+(multMI-1)]+((lset-setBound[1])/nAttParAxis)*nDegFreedomAtt+(localAstroMax-offsetAttParam);
                            }
                            else
                                j++;
                        }
                    sum=sum+systemMatrix[l]*nullSpaceVect[j];
                    
                    double NSVal;
                    if(lset<setBound[1]){
                            NSVal=(matrixIndex[multMI*i]+lset)+ic*nStar*nAstroPSolved;
                            sumVer=sumVer+1.0*NSVal;
                    }
                    else{
                                if(ic==0 || ic==3){
                                
                                if(lset<setBound[1]+4)
                                    NSVal=(matrixIndex[multMI*i+1]-nStar*nAstroPSolved)+lset-setBound[1]+ic*nDegFreedomAtt;
                                else
                                    NSVal=0;
                                }
                               
                            if(ic==1 || ic==4){
                                if(lset>=setBound[1]+4 && lset<setBound[1]+8)
                                    NSVal=(matrixIndex[multMI*i+1]-nStar*nAstroPSolved)+ic*nDegFreedomAtt+lset-(setBound[1]+4);
                                else
                                    NSVal=0;
                            }
                            if(ic==2 || ic==5){
                                if(lset>=setBound[1]+8)
                                    NSVal=(matrixIndex[multMI*i+1]-nStar*nAstroPSolved)+ic*nDegFreedomAtt+lset-(setBound[1]+8);
                                else
                                    NSVal=0;
                            }
                            sumVer=sumVer+1.0*NSVal;
                    }
                     lset++;
                }//for(l
/*                    if(i>=nElements){
                            printf("ERROR: PE=%d i=%d   greather than nElements=%ld\n",myid,i ,nElements);
                            continue;
                    }*/
               if(sumVer != chkSumVer){
                   printf("ERROR: PE=%d NullSapce Equation ic=%d, sumVer[%d]=%d and chkSumVer=%ld are not equal\n",myid,ic,i,sumVer,    chkSumVer);
                   MPI_Abort(MPI_COMM_WORLD, 0);
                   exit(1);
               }
                productNS[i]=sum;
        }//for i
        ////////////   in case of ExtConstr
/*  TO BE DECIDED orthogonalty of nullspace base vectors
         for(i2=0;i2<nEqExtConstr;i2++ ){
         sum=0.0;
         for(j2=0;j2<numOfExtStar*nAstroPSolved;j2++){
         sum+=(systemMatrix[mapNcoeff[myid]+j2+i2*nOfElextObs]/extConstrW)*nullSpaceVect[j2];
         }
         for(na=0;na<nAttAxes;na++)
         for(j2=0;j2<numOfExtAttCol;j2++){
         sum+=(systemMatrix[mapNcoeff[myid]+numOfExtStar*nAstroPSolved+j2+i2*nOfElextObs+na*numOfExtAttCol]/extConstrW)*nullSpaceVect[localAstroMax+startingAttColExtConstr+na*nDegFreedomAtt+j2];
         }
         productNS[mapNoss[myid]+i2]+=sum;
         }//for i2
         
         double *kcopy;
         kcopy=(double *) calloc(nEqExtConstr, sizeof(double));
         MPI_Allreduce(&productNS[mapNoss[myid]],kcopy,nEqExtConstr,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
         for(int i3=0;i3<nEqExtConstr;i3++)
         {
         productNS[mapNoss[myid]+i3]=kcopy[i3];
         }
         free(kcopy);
         ////////////   in case of ExtConstr
         */
}//pragma
        double normLoc;
        normLoc=cblas_dnrm2(mapNoss[myid],productNS,1);
        double normLoc2=normLoc*normLoc;
        double nrmGlob;
        MPI_Allreduce(&normLoc2, &nrmGlob,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
        results.vectNorm[ic]=sqrt(nrmGlob);
        double localMin=productNS[0];
        double localMax=productNS[0];
        double globalSum=0.0;
        double localSum=0.0;
        for(int j=0;j<comlsqr.mapNoss[myid];j++){
//            if(myid==39 && j==1456748) printf("TP0.1 PE=%d   localSum=%lf localMin=%lf localMax=%lf productNS[1456748]=%lf\n",myid,localSum,localMin,localMax,productNS[1456748]);
           if(localMin>productNS[j]) localMin=productNS[j];
            if(localMax<productNS[j]) localMax=productNS[j];
            localSum+=productNS[j];
//                if(myid==39 && j==1456748) printf("TP0.2 PE=%d   localSum=%lf localMin=%lf localMax=%lf productNS[1456748]=%lf\n",myid,localSum,localMin,localMax,productNS[1456748]);
        }
//       if(myid==39) printf("TP1 PE=%d   localSum=%lf localMin=%lf localMax=%lf productNS[1456748]=%lf\n",myid,localSum,localMin,localMax,productNS[1456748]);
       MPI_Allreduce(&localMin,&results.compMin[ic], 1, MPI_DOUBLE, MPI_MIN,MPI_COMM_WORLD );
       MPI_Allreduce(&localMax,&results.compMax[ic], 1, MPI_DOUBLE, MPI_MAX,MPI_COMM_WORLD );
       MPI_Allreduce(&localSum,&globalSum, 1, MPI_DOUBLE, MPI_SUM,MPI_COMM_WORLD );
       double avg=globalSum/comlsqr.nobs;
       results.compAvg[ic]=avg;

       double localsqrsum=0;
       double globalsqrsum=0;
       for(int j=0;j<comlsqr.mapNoss[myid];j++)
               localsqrsum+=productNS[j]*productNS[j];
               
        MPI_Allreduce(&localsqrsum,&globalsqrsum, 1, MPI_DOUBLE, MPI_SUM,MPI_COMM_WORLD );
        results.compVar[ic]=sqrt((globalsqrsum/comlsqr.nobs)-avg*avg);
 
        if(myid==0){
               FILE *fpNS;
               char icstr[10];
               char *fileName;
               fileName=(char *) malloc(512);
               sprintf(icstr,"%1d",ic);
               strcpy(fileName, "Gsr_NullSpaceProduct_");
               strcat(fileName, icstr);
               strcat(fileName, ".dat");
               fpNS=fopen(fileName,"w");
               if (!fpNS)
               {
                   printf("PE=%d Error while open %s\n",myid,fileName);
                   MPI_Abort(MPI_COMM_WORLD,1);
                   exit(EXIT_FAILURE);
               }
               long t1=comlsqr.nobs/100;
               int bins=MIN(1000,t1);
               int *isto;
               double binWidth=(results.compMax[ic]-results.compMin[ic])/bins;
                 int ix;
                isto=(int *) calloc(bins, sizeof(int));
                 for( int j=0;j<mapNoss[myid];j++){
                     if(productNS[j]==results.compMax[ic]){
                         ix=bins-1;
                         isto[ix]++;
                         continue;
                     }
                     ix=(productNS[j]-results.compMin[ic])/binWidth;
                     if(ix >=bins){
                         printf("Error while computing local myid=0 ix=%d on cknullSpace at element %d on equation %d. prodNS=%lf greater than max=%lf\n",ix,j,ic,productNS[j],results.compMax[ic]);
                         MPI_Abort(MPI_COMM_WORLD,1);
                         exit(EXIT_FAILURE);
                        
                     }
                     isto[ix]++;
                 }
              for(int k=1;k<nproc;k++){
                   MPI_Recv(productNS, mapNoss[k], MPI_DOUBLE, k, 0, MPI_COMM_WORLD,&status);
//                  if(k==39) printf("TP2 from PE=%k  productNS[1456748]=%lf\n",k,productNS[1456748]);
                  for( int j=0;j<mapNoss[k];j++){
                      if(productNS[j]==results.compMax[ic]){
                          ix=bins-1;
                          isto[ix]++;
                          continue;
                      }
                      ix=(productNS[j]-results.compMin[ic])/binWidth;
                      if(ix >=bins){
                          printf("Error while computing remote peid=%d ix=%d on cknullSpace at element %d on equation %d. prodNS=%lf greater than max=%lf\n",k,ix,j,ic,productNS[j],results.compMax[ic]);
                          MPI_Abort(MPI_COMM_WORLD,1);
                          exit(EXIT_FAILURE);
                          
                      }
                      isto[ix]++;
                  }

               }
              fprintf(fpNS, "%le\n%le\n%le\n%le\n%le\n%d\n%le\n",results.compMin[ic],results.compMax[ic],results.compAvg[ic],results.compVar[ic],
                      results.vectNorm[ic],bins,binWidth);
            for(int j=0;j<bins;j++)
                fprintf(fpNS, "%d\n",isto[j]);
               fclose(fpNS);
               free(fileName);
                free(isto);
                seconds[1] = time(NULL);
                tot_sec = seconds[1] - seconds[0];
                 printf("Istogram nullSpace file for Nullspace row %d, in %ld sec\n",ic,tot_sec);

             } else{
//                 if(myid==39) printf("TP3 PE=%d  productNS[1456748]=%lf\n",myid,productNS[1456748]);

               MPI_Send(productNS, mapNoss[myid], MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);
        }//if myid==0
        
        } //for ic
    free(nullSpaceFPN);
    free(productNS);
    free(prodNS);
    free(nullSpaceVect);
    return results;
}
char* subString (const char* input, int offset, int len, char* dest)
{
    int input_len = strlen (input);
    
    if (offset + len > input_len)
    {
        return NULL;
    }
    
    strncpy (dest, input + offset, len);
    return dest;
}

/*--------------------------------------------------------------------------*/
// Compute the value of the shifted Legendre polinomial of degree deg. Here we
// just need deg<3.
double legendre(int deg, double x) {
	double res;
	
	switch(deg) {
	case 0:
		res = 1.0;
		break;
	case 1:
		res = 2.0*(x-0.5);
		break;
	case 2:
		res = 6.0*(x-0.5)*(x-0.5)-0.5;
		break;
	default:
		res = -1.0;
		break;
}
	
	return res;
}

//void computeInstrContr (int lsInstrFlag, int ssInstrFlag, long nFoVs,long nCCDs, long nPixelColumns,long nTimeIntervals,double * instrCoeffConstr,int * instrColsConstr,int * instrConstrIlung)
int computeInstrConstr (struct comData comlsqr,double * instrCoeffConstr,int * instrColsConstr,int * instrConstrIlung)
{
    ////////////////// Writing instrConstrRows_xxxx.bin file
// There are 1 AL + 2 AC constraint equations for each time interval for the large scale parameters (total=3*nTimeIntervals)
// There are 1 AL + 1 AC constraint equations for each CCD and Legendre polinomial degree (total=2*nCCDs*3)
// The equations are described by three arrays: instrCoeffConstr, instrColsConstr, instrConstrIlung, which contain the
// coefficients, the column indexes and the length of the non-zero coefficients of each equation respectively.
// MUST BE ALLOCATED THE FOLLOWING VECTORS:
//instrCoeffConstr=(double *) calloc(nElemIC, sizeof(double));
//instrColsConstr=(int *) calloc(nElemIC, sizeof(int));
//instrConstrIlung=(int *) calloc(nOfInstrConstr, sizeof(int));

    int lsInstrFlag=comlsqr.lsInstrFlag;
    int ssInstrFlag=comlsqr.ssInstrFlag;
    long nFoVs=1+comlsqr.instrConst[0];
    long nCCDs=comlsqr.instrConst[1];
    long nPixelColumns=comlsqr.instrConst[2];
    long nTimeIntervals=comlsqr.instrConst[3];
    long cCDLSAACZP=comlsqr.cCDLSAACZP;
    
int nElemICLSAL =comlsqr.nElemICLSAL;
int nElemICLSAC =comlsqr.nElemICLSAC;
int nElemICSS = comlsqr.nElemICSS;
int nOfInstrConstr = comlsqr.nOfInstrConstr;
int nElemIC = comlsqr.nElemIC;
int counterElem=0;
int counterEqs=0;
int elemAcc=0;

if(lsInstrFlag){
    // generate large scale constraint eq. AL
    for(int i=0; i<nTimeIntervals; i++) {
        instrConstrIlung[counterEqs] = nElemICLSAL;
        elemAcc+=nElemICLSAL;
        counterEqs++;
        for(int j=0; j<nFoVs; j++) {
            for(int k=0; k<nCCDs; k++) {
                instrCoeffConstr[counterElem] = 1.0;
                instrColsConstr[counterElem] = comlsqr.offsetCdelta_eta + j*nCCDs*nTimeIntervals+k*nTimeIntervals+i;
                counterElem++;
            }
        }
    }
    // generate large scale constraint eq. AC
    for(int i=0; i<nTimeIntervals; i++) {
        for(int j=0; j<nFoVs; j++) {
            instrConstrIlung[counterEqs] = nElemICLSAC;
            elemAcc+=nElemICLSAC;
            counterEqs++;
            for(int k=0; k<nCCDs; k++) {
                instrCoeffConstr[counterElem] = 1.0;
                instrColsConstr[counterElem] = comlsqr.offsetCdelta_zeta + j*nCCDs*nTimeIntervals+k*nTimeIntervals+i;
                counterElem++;
            }
        }
    }