/*--------------------------------------------*/
/*                                            */
/* (c) Sebastian Hiller, 2004-2007            */
/*                                            */
/*--------------------------------------------*/

#include <stdio.h>
#include <math.h>
#include <string.h>

/*--------------------------------------------*/
/* Before compilation, define these variables:*/
/*                                            */

#define LOCATIONGAPRO "/home/nmrsu/gapro_098_linux"

#define PLATFORM "XWINNMR 3.5"    /* or "TOPSPIN 2.0", "TOPSPIN 2.1"*/

#define TEXTEDITOR "/home/nmrsu/nedit"

/*                                            */
/*                                            */
/*--------------------------------------------*/
/*  version 01 SH 04dec07                     */
/*  version 02 SH 23jan08 fullrun enabled     */
/*  version 03 SH 24sep08 new format          */
/*--------------------------------------------*/

#define MAXPROJECTIONS 500
#define STRING 256
#define MAXDIMENSIONS 10

/*----------------------definition of variables--------------------*/

FILE *spectrafile, *spectrafile2, *prosafile, *anglefile, *logfile, *communfile;
FILE *megaprosa, *errorfile, *controlfile, *pprogfile;

char mainpath[STRING], greekanglenames[25][10], gapro[STRING], spectra[STRING], spectra2[STRING], spectrum[STRING], expfilename[STRING];
char prosa[STRING], communication[STRING], log[STRING], execute[500], anglefname[200], errorfname[STRING], spectrum2[STRING];
char editor[STRING], parfilename[STRING], pulprog[100], pprogname[STRING], experiment[30];
char freqs[3][10], carriers[MAXDIMENSIONS][10], cnstnum[6], ttemp[20], megaprosname[STRING], zgadd[5], zgadd2[2];
char dimusedlett[MAXPROJECTIONS][MAXDIMENSIONS], zgoptns[STRING], phaseoverlap[30], phasesused[30];
char dimusedlett3D[MAXPROJECTIONS][3][MAXDIMENSIONS], errorstring[256], inputline[STRING], byteorder[STRING];

int validanswer, convergence, fixednumproj, td2, usefile, trosydim, dimensionality;
int apsycomplete, numofangles[3], numofcombidim[3], totalnumproj, td1, trosy, pts2complex, pts1complex;
int numofp, transient, decimm, basexpno, numberscans, number, dimproj, dataset, pts3complex, td3;
int i, j, k, kk, countter, count2, fullrun, setuponly, finalsi2, limit, poexists;
int si1, si2, si3, stsi2, stsr2, report, numofprojsofar, finished, peaksfound[MAXPROJECTIONS*4];
int ptomeasure[MAXPROJECTIONS], nondangles[MAXPROJECTIONS], plusminus[MAXPROJECTIONS][MAXDIMENSIONS];
int dimused[MAXPROJECTIONS][MAXDIMENSIONS],parmode,combidim[MAXPROJECTIONS][MAXDIMENSIONS][3];
int nondangles3D[MAXPROJECTIONS][3], plusminus3D[MAXPROJECTIONS][MAXDIMENSIONS][3];
int dimused3D[MAXPROJECTIONS][MAXDIMENSIONS][3];

float angles[MAXPROJECTIONS][MAXDIMENSIONS], angles3D[MAXPROJECTIONS][MAXDIMENSIONS][3],SWhertz, sweept, sweepu, sweepv;
float sweepw, sweepx, sweepy, temp, sweep1, sweep3, pphhase, pphhase2, sweep_a[MAXDIMENSIONS];

double groupdelay, freq[3], carrier[MAXDIMENSIONS];

const int addtheorem[8][8]={
{1, 1, 1, 1, 1, 1, 1, 1},
{1,-1, 1,-1, 1,-1, 1,-1},
{1, 1,-1,-1, 1, 1,-1,-1},
{1,-1,-1, 1, 1,-1,-1, 1},
{1, 1, 1, 1,-1,-1,-1,-1},
{1,-1, 1,-1,-1, 1,-1, 1},
{1, 1,-1,-1,-1,-1, 1, 1},
{1,-1,-1, 1,-1, 1, 1,-1}
};

const int minuss[8][4]={
{0,0,0,0},
{0,1,0,0},
{0,0,1,0},
{0,1,1,0},
{0,0,0,1},
{0,1,0,1},
{0,0,1,1},
{0,1,1,1}
};

const int tppi[32][6]={
{0,0,0,0,0,0},
{0,0,0,0,1,1},
{0,0,0,1,0,1},
{0,0,0,1,1,0},
{0,0,1,0,0,1},
{0,0,1,0,1,0},
{0,0,1,1,0,0},
{0,0,1,1,1,1},
{0,1,0,0,0,1},
{0,1,0,0,1,0},
{0,1,0,1,0,0},
{0,1,0,1,1,1},
{0,1,1,0,0,0},
{0,1,1,0,1,1},
{0,1,1,1,0,1},
{0,1,1,1,1,0},
{1,0,0,0,0,1},
{1,0,0,0,1,0},
{1,0,0,1,0,0},
{1,0,0,1,1,1},
{1,0,1,0,0,0},
{1,0,1,0,1,1},
{1,0,1,1,0,1},
{1,0,1,1,1,0},
{1,1,0,0,0,0},
{1,1,0,0,1,1},
{1,1,0,1,0,1},
{1,1,0,1,1,0},
{1,1,1,0,0,1},
{1,1,1,0,1,0},
{1,1,1,1,0,0},
{1,1,1,1,1,1}
};

const int receiver[32]={
0,2,2,2,2,2,2,0,2,2,2,0,2,0,0,0,
2,2,2,0,2,0,0,0,2,0,0,0,0,0,0,2
};

/*----------------------START of program------------------------------*/

/*----------------------Initialization--------------------------------*/

sprintf(gapro,LOCATIONGAPRO);  /*location of gapro*/

strcpy(greekanglenames[1],"alpha");
strcpy(greekanglenames[2],"beta");
strcpy(greekanglenames[3],"gamma");
strcpy(greekanglenames[4],"delta");
strcpy(greekanglenames[5],"epsilon");
strcpy(greekanglenames[6],"zeta");
strcpy(greekanglenames[7],"eta");
strcpy(greekanglenames[8],"theta");
strcpy(greekanglenames[9],"iota");
strcpy(greekanglenames[10],"kappa");
strcpy(greekanglenames[11],"lambda");
strcpy(greekanglenames[12],"my");
strcpy(greekanglenames[13],"ny");
strcpy(greekanglenames[14],"xi");
strcpy(greekanglenames[15],"omikron");
strcpy(greekanglenames[16],"pi");
strcpy(greekanglenames[17],"rho");
strcpy(greekanglenames[18],"sigma");
strcpy(greekanglenames[19],"tau");
strcpy(greekanglenames[20],"ypsilon");
strcpy(greekanglenames[21],"phi");
strcpy(greekanglenames[22],"chi");
strcpy(greekanglenames[23],"psi");
strcpy(greekanglenames[24],"omega");

zgadd[0]=' ';
zgadd[1]='-';
zgadd[2]='D';
zgadd[3]=' ';
zgadd[4]='\0';

zgadd2[0]=' ';
zgadd2[1]='\0';

validanswer=0;
apsycomplete=0;
usefile=0;
convergence=0;
fixednumproj=0;
totalnumproj=0;
trosy=0;
fullrun=0;
setuponly=0;
numofprojsofar=0;
transient=0;
poexists=0;
report=1;
parmode=0;

GETCURDATA
FETCHPAR("PARMODE",&parmode);
FETCHPAR("PULPROG",pulprog);
FETCHPAR("TD",&td2)
pts2complex =td2/2;
FETCHPAR1("TD",&td1)
pts1complex=td1/2;

if(strcmp(PLATFORM,"XWINNMR 3.5")==0){
  strcpy(byteorder,"swap");
}
else if(strcmp(PLATFORM,"TOPSPIN 2.0")==0){
  strcpy(byteorder,"integer");
}
else if(strcmp(PLATFORM,"TOPSPIN 2.1")==0){
  strcpy(byteorder,"integer");
}
else {
  STOPMSG("Platform not defined")
}


if(strcmp(PLATFORM,"XWINNMR 3.5")==0){
  FETCHPAR("DECIM",&decimm)
  /*----------following values are only valid for dspfvs==12------*/
  if(decimm==2) groupdelay=46.311;
  if(decimm==3) groupdelay=36.530;
  if(decimm==4) groupdelay=47.870;
  if(decimm==6) groupdelay=50.229;
  if(decimm==8) groupdelay=53.289;
  if(decimm==12) groupdelay=69.551;
  if(decimm==16) groupdelay=71.600;
  if(decimm==24) groupdelay=70.167;
  if(decimm==32) groupdelay=72.138;
  if(decimm==48) groupdelay=70.528;
  if(decimm==64) groupdelay=72.348;
  if(decimm==96) groupdelay=70.700;
  if(decimm==128) groupdelay=72.524;
}
else if(strcmp(PLATFORM,"TOPSPIN 2.0")==0){
  FETCHPAR("GRPDLY",&groupdelay)
}
else if(strcmp(PLATFORM,"TOPSPIN 2.1")==0){
  FETCHPAR("GRPDLY",&groupdelay)
}
else {
  STOPMSG("Platform not defined")
}
pphhase2=groupdelay*360.0;

FETCHPAR("SWH",&SWhertz)
FETCHPAR("PHC0",&pphhase)
FETCHPAR("NS",&numberscans)
FETCHPAR1("SI",&si1)
FETCHPAR("SI",&si2)
FETCHPAR("STSI",&stsi2)
FETCHPAR("STSR",&stsr2)
if(stsi2==0) {
  finalsi2=si2;
}
else {
  finalsi2=stsi2;
}
/*
editor=getenv("EDITOR");
if (editor==NULL) {
  STOPMSG("please define environment variable: EDITOR")
}
*/
strcpy(editor,TEXTEDITOR);

sprintf(mainpath,"%s/data/%s/nmr/%s/%d",disk,user,name,expno);
sprintf(spectra,"%s/spectra.gap",mainpath);
sprintf(spectra2,"%s/spectra.gap_prosa",mainpath);
sprintf(log,"%s/APSY.log",mainpath);
sprintf(errorfname,"%s/error.gap",mainpath);
sprintf(communication,"%s/term.gap",mainpath);
sprintf(anglefname,"%s/angles.dat",mainpath);
sprintf(expfilename,"%s/experiment.gap",mainpath);
sprintf(parfilename,"%s/parameter.gap",mainpath);
sprintf(pprogname,"%s/exp/stan/nmr/lists/pp/%s",disk,pulprog);
sprintf(megaprosname,"%s/all_apsy.pro",mainpath);

if((pprogfile=fopen(pprogname,"r"))==NULL) {
  STOPMSG("can not open pulse program\n")
}
fscanf(pprogfile,";;Dimensionality:");
fscanf(pprogfile,"%d",&dimensionality);
fscanf(pprogfile,"\n");
fscanf(pprogfile,";;Experiment:");
fscanf(pprogfile,"%s",experiment);
fscanf(pprogfile,"\n");
fscanf(pprogfile,";;Phaseoverlap:");
fscanf(pprogfile,"%s",phaseoverlap);
fscanf(pprogfile,"\n");
fscanf(pprogfile,";;trosy:");
fscanf(pprogfile,"%d",&trosydim);
fscanf(pprogfile,"\n");
fscanf(pprogfile,";;DimProjections:");
fscanf(pprogfile,"%d",&dimproj);
fscanf(pprogfile,"\n");
if(dimproj!=parmode+1){
  STOPMSG("Inconsistency between Parmode and DimProj\n")
}
fscanf(pprogfile,";;Hfreq:");
fscanf(pprogfile,"%s",freqs[0]);
fscanf(pprogfile,"\n");
fscanf(pprogfile,";;Cfreq:");
fscanf(pprogfile,"%s",freqs[1]);
fscanf(pprogfile,"\n");
fscanf(pprogfile,";;Nfreq:");
fscanf(pprogfile,"%s",freqs[2]);
for(i=1;i<=dimensionality;i++){
  fscanf(pprogfile,"\n");
  fscanf(pprogfile,";;Carrier");
  fscanf(pprogfile,"%d",&number);
  if(number!=i){
    STOPMSG("wrong header format");
  }
  fscanf(pprogfile,":");
  fscanf(pprogfile,"%s",carriers[i]);
  if(strcmp(carriers[i],"CNST")==0) {
    fscanf(pprogfile,"%s",cnstnum);
    sprintf(ttemp,"%s %s",carriers[i],cnstnum);
    sprintf(carriers[i],"%s",ttemp);
  }
}
fclose(pprogfile);

if((dimensionality<3)||(dimensionality>7)){
  STOPMSG("This dimensionality is not supported");
}
if((dimproj<2)||(dimproj>3)){
  STOPMSG("This dimensionality of projections is not supported");
}

if(dimproj==3){
  FETCHPAR3("TD",&td3)
  pts3complex=td3/2;
  FETCHPAR3("SI",&si3)
}

if((megaprosa=fopen(megaprosname,"w"))==NULL) {
  STOPMSG("can not open multi prosa macro\n")
}
fprintf(megaprosa,"linpred  := n\n");
fprintf(megaprosa,"baseline2  := y\n");
fprintf(megaprosa,"baseline3  := y\n");
fprintf(megaprosa,"mirrlp   := n\n");
fprintf(megaprosa,"wfunction   := \"sin 75\"\n");
fprintf(megaprosa,"shapes     = 20\n\n");
fprintf(megaprosa,"ptspredict = 32\n\n");
fprintf(megaprosa,"phcor    = %4.1f\n\n",pphhase);
fprintf(megaprosa,"ptindi = %d\n\n", si1);
fprintf(megaprosa,"ptindi3 = %d\n\n", si3);
fprintf(megaprosa,"folder := \"choose_a_name\"\n\n");
fprintf(megaprosa,"system mkdir /nmr/data/%s/%s/$folder\n",user,name);
fprintf(megaprosa,"system cp experiment.gap /nmr/data/%s/%s/$folder/experiment.gap\n",user,name);
fprintf(megaprosa,"system cp parameter.gap /nmr/data/%s/%s/$folder/parameter.gap\n",user,name);
fprintf(megaprosa,"system cp spectra.gap_prosa /nmr/data/%s/%s/$folder/spectra.gap\n\n",user,name);
fclose(megaprosa);

if(dimproj==2) {
  numofangles[1]=dimensionality-2;
  numofcombidim[1]=dimensionality-1;
}
else {
  if(dimensionality==4){
    numofangles[1]=1;
    numofangles[2]=0;
  }
  else if(dimensionality==5){
    numofangles[1]=1;
    numofangles[2]=1;
  }
  else {
    STOPMSG("not supported code 1");
  }
  numofcombidim[1]=numofangles[1]+1;
  numofcombidim[2]=numofangles[2]+1;
}

for(i=0;i<dimensionality;i++){
  if(phaseoverlap[i]!='0') poexists++;
}

for(i=0;i<3;i++){
  if(strncmp(freqs[i],"BF",2)==0){
    FETCHPAR(freqs[i], freq+i)
  }
  else {
    STOPMSG("This freq string is not supported");
  }
}

for(i=1;i<=dimensionality;i++){
  if(strncmp(carriers[i],"O",1)==0){
    FETCHPAR(carriers[i], carrier+i)
  }
  else if(strncmp(carriers[i],"CNST",4)==0){
    FETCHPAR(carriers[i], &temp)
    carrier[i]=((double) temp);
  }
  else {
    STOPMSG("This carrier string is not supported");
  }
}

if((controlfile=fopen(expfilename, "w"))==NULL) {
  STOPMSG("can not create experiment.gap.\n")
}
fprintf(controlfile,"DimensionalityProjections:\t%d\n",dimproj);
fprintf(controlfile,"DimensionalitySpace:\t\t%d\n",dimensionality);
for(i=1;i<=dimensionality;i++){
  fprintf(controlfile,"Carrier%d:\t\t\t%f\n",i,carrier[i]);
}
fprintf(controlfile,"Hfreq:\t\t\t\t%f\n",freq[0]);
fprintf(controlfile,"Nfreq:\t\t\t\t%f\n",freq[2]);
fprintf(controlfile,"Cfreq:\t\t\t\t%f\n",freq[1]);
fprintf(controlfile,"Experiment:\t\t\t%s\n",experiment);
fclose(controlfile);

if((spectrafile=fopen(spectra,"w"))==NULL) {
  STOPMSG("can not open spectra file\n")
}
fclose(spectrafile);

if((spectrafile2=fopen(spectra2,"w"))==NULL) {
  STOPMSG("can not open spectra file 2\n")
}
fclose(spectrafile2);

if((logfile=fopen(log,"w"))==NULL){
  STOPMSG("can not open logfile\n");
}
fclose(logfile);

sprintf(execute,"%s %s &",editor,log);
system(execute);

if((communfile=fopen(communication,"r"))!=NULL){
  fclose(communfile);
  sprintf(execute,"rm &s",communication);
  system(execute);
  ssleep(0.1);
}

while(validanswer==0){
  if(poexists){
    GETINT("APSYpo. Choose an option:\n Full run (1)\n Set up experiments (2)\n Cancel (0)",i1)
  }
  else {
    GETINT("APSY. Choose an option:\n Full run (1)\n Set up experiments (2)\n Cancel (0)",i1)
  }
  if(i1==1){
    validanswer=1;
    fullrun=1;
  }
  else if(i1==2){
    validanswer=1;
    setuponly=1;
  }
  else if(i1==0){
    STOP
  }
}

setuponly=1;
validanswer=0;
i1=0;

if(fullrun){
  if((controlfile=fopen(parfilename, "r"))==NULL) {
    sprintf(execute,"%s -CreateParamTemplate %s %d", gapro, mainpath, dimensionality);
    system(execute);
    finished=0;
    while(finished==0){
      ssleep(0.1);
      if((communfile=fopen(communication,"r"))!=NULL){
        finished=1;
        fclose(communfile);
      }
      if((errorfile=fopen(errorfname,"r"))!=NULL){
        fclose(errorfile);
        sprintf(execute,"rm %s",errorfname);
        system(execute);
        ssleep(0.1);
        STOPMSG("Unusual termination of GAPRO");
      }
    }
    sprintf(execute,"rm %s",communication);
    system(execute);
    ssleep(0.1);
    sprintf(execute,"%s %s &",editor,parfilename);
    system(execute);
    while(validanswer==0){
      GETINT("A template parameter.gap has been created.\nPlease specify parameters and \n continue (1)\n cancel (2)",i1)
      if(i1==1){
        validanswer=1;
      }
      else if(i1==2){
        STOP
      }
    }
  }
  else {
    sprintf(execute,"%s %s &",editor,parfilename);
    system(execute);
    while(validanswer==0){
      GETINT("The file parameter.gap is used.\nPlease check parameters and \n continue (1)\n cancel (2)",i1)
      if(i1==1){
        validanswer=1;
      }
      else if(i1==2){
        STOP
      }
    }
  }
}

validanswer=0;
i1=0;

while(validanswer==0){
  GETINT("Choose an option:\n Use a list of angles (1)\n Enter one set of angles manually (2)",i1)
  if(i1==1){
    validanswer=1;
    usefile=1;
  }
  else if(i1==2){
    validanswer=1;
    fixednumproj=1;
    totalnumproj=1;
  }
}

/*-------------------------------manually entering of angles------------------------------*/
if(usefile==0){
  for(i=1;i<=numofangles[1];i++) {
    GETFLOAT(greekanglenames[i],angles[0][i-1])
  }
}

/*----------------------------reading of anglefile---------------------------------------*/
if(usefile){
  if((anglefile=fopen(anglefname, "r"))==NULL){
    if(dimproj==2){
      sprintf(execute,"%s -XWINMatchingAngles %s %d", gapro, mainpath, dimensionality);
      system(execute);
      finished=0;
      while(finished==0){
        ssleep(0.1);
        if((communfile=fopen(communication,"r"))!=NULL){
          finished=1;
          fclose(communfile);
        }
        if((errorfile=fopen(errorfname,"r"))!=NULL){
          fclose(errorfile);
          sprintf(execute,"rm %s",errorfname);
          system(execute);
          ssleep(0.1);
          STOPMSG("Unusual termination of GAPRO");
        }
      }
      sprintf(execute,"rm %s",communication);
      system(execute);
      ssleep(0.1);
    }
    else{
    if((anglefile=fopen(anglefname, "w"))==NULL) {
      STOPMSG("Can not create angle.dat file. Please check your permissions\n")
    }
      if(dimensionality==4){
        fprintf(anglefile,"dim1=12 alpha = 0.0 dim2=3\n");
        fprintf(anglefile,"dim1=13 alpha = 0.0 dim2=2\n");
        fprintf(anglefile,"dim1=23 alpha = 0.0 dim2=1\n");
        fprintf(anglefile,"dim1=12 alpha = 30.0 dim2=3\n");
        fprintf(anglefile,"dim1=13 alpha = 40.0 dim2=2\n");
        fprintf(anglefile,"dim1=23 alpha = 50.0 dim2=1\n");
      }
      else if(dimensionality==5){
        fprintf(anglefile,"dim1=12 alpha =  0.0 dim2=34  alpha =  0.0\n");
        fprintf(anglefile,"dim1=12 alpha =  0.0 dim2=34  alpha = 90.0\n");
        fprintf(anglefile,"dim1=12 alpha =  0.0 dim2=34  alpha = 45.0\n");
        fprintf(anglefile,"dim1=12 alpha = 90.0 dim2=34  alpha =  0.0\n");
        fprintf(anglefile,"dim1=12 alpha = 90.0 dim2=34  alpha = 90.0\n");
        fprintf(anglefile,"dim1=12 alpha = 90.0 dim2=34  alpha = 45.0\n");
        fprintf(anglefile,"dim1=12 alpha = 45.0 dim2=34  alpha =  0.0\n");
        fprintf(anglefile,"dim1=12 alpha = 45.0 dim2=34  alpha = 90.0\n");
        fprintf(anglefile,"dim1=12 alpha = 45.0 dim2=34  alpha = 45.0\n");
        fprintf(anglefile,"dim1=13 alpha =  0.0 dim2=24  alpha =  0.0\n");
        fprintf(anglefile,"dim1=13 alpha =  0.0 dim2=24  alpha = 90.0\n");
        fprintf(anglefile,"dim1=13 alpha =  0.0 dim2=24  alpha = 45.0\n");
        fprintf(anglefile,"dim1=13 alpha = 90.0 dim2=24  alpha =  0.0\n");
        fprintf(anglefile,"dim1=13 alpha = 90.0 dim2=24  alpha = 90.0\n");
        fprintf(anglefile,"dim1=13 alpha = 90.0 dim2=24  alpha = 45.0\n");
        fprintf(anglefile,"dim1=13 alpha = 45.0 dim2=24  alpha =  0.0\n");
        fprintf(anglefile,"dim1=13 alpha = 45.0 dim2=24  alpha = 90.0\n");
        fprintf(anglefile,"dim1=13 alpha = 45.0 dim2=24  alpha = 45.0\n");
        fprintf(anglefile,"dim1=14 alpha =  0.0 dim2=23  alpha =  0.0\n");
        fprintf(anglefile,"dim1=14 alpha =  0.0 dim2=23  alpha = 90.0\n");
        fprintf(anglefile,"dim1=14 alpha =  0.0 dim2=23  alpha = 45.0\n");
        fprintf(anglefile,"dim1=14 alpha = 90.0 dim2=23  alpha =  0.0\n");
        fprintf(anglefile,"dim1=14 alpha = 90.0 dim2=23  alpha = 90.0\n");
        fprintf(anglefile,"dim1=14 alpha = 90.0 dim2=23  alpha = 45.0\n");
        fprintf(anglefile,"dim1=14 alpha = 45.0 dim2=23  alpha =  0.0\n");
        fprintf(anglefile,"dim1=14 alpha = 45.0 dim2=23  alpha = 90.0\n");
        fprintf(anglefile,"dim1=14 alpha = 45.0 dim2=23  alpha = 45.0\n");
       }
      else {
        STOPMSG("not supported code 4");
      }
      fclose(anglefile);
    }
    sprintf(execute,"%s %s &",editor,anglefname);
    system(execute);
    validanswer=0;
    i1=0;
    while(validanswer==0){
      GETINT("A template angles.dat has been created.\nPlease specify angles and \n continue (1)\n cancel (2)",i1)
      if(i1==1){
        validanswer=1;
      }
      else if(i1==2){
        STOP
      }
    }
  }
  else {
    fclose(anglefile);
    sprintf(execute,"%s %s &",editor,anglefname);
    system(execute);
    validanswer=0;
    i1=0;
    while(validanswer==0){
      GETINT("The file angles.dat is used.\nPlease check the angles and \n continue (1)\n cancel (2)",i1)
      if(i1==1){
        validanswer=1;
      }
      else if(i1==2){
        STOP
      }
    }
  }
  if((anglefile=fopen(anglefname, "r"))==NULL){
    STOPMSG("can not read anglefile\n")
  }
  i=0;
  if(dimproj==2){
    while(fgets(inputline,STRING,anglefile)!=NULL){
      if(strncmp(inputline,"#",1)!=0){
        if(numofangles[1]==1){
          sscanf(inputline,"%f",angles[i]+0);
	}
        else if(numofangles[1]==2){
          sscanf(inputline,"%f %f",angles[i]+0,angles[i]+1);
	}
        else if(numofangles[1]==3){
          sscanf(inputline,"%f %f %f",angles[i]+0,angles[i]+1,angles[i]+2);
	}
        else if(numofangles[1]==4){
          sscanf(inputline,"%f %f %f %f",angles[i]+0,angles[i]+1,angles[i]+2,angles[i]+3);
	}
        else if(numofangles[1]==5){
          sscanf(inputline,"%f %f %f %f %f",angles[i]+0,angles[i]+1,angles[i]+2,angles[i]+3,angles[i]+4);
	}
        else if(numofangles[1]==6){
          sscanf(inputline,"%f %f %f %f %f %f",angles[i]+0,angles[i]+1,angles[i]+2,angles[i]+3,angles[i]+4,angles[i]+5);
	}
	else {
          STOPMSG("not implemented code 38");
	}
        i++;
      }
    }
  }
  else if(dimproj==3){
    if(dimensionality==4){
      while(fscanf(anglefile,"dim1=%1d%1d",&(combidim[i][0][1]),&(combidim[i][1][1]))==2){
        fscanf(anglefile," ");
        fscanf(anglefile,"alpha");
        fscanf(anglefile," = ");
        fscanf(anglefile,"%f",angles3D[i][0]+1);
        fscanf(anglefile,"   ");
        fscanf(anglefile,"dim2=%1d",&(combidim[i][0][2]));
        fscanf(anglefile,"\n");
        i++;
      }
    }
    else if(dimensionality==5){
      while(fscanf(anglefile,"dim1=%1d%1d",&(combidim[i][0][1]),&(combidim[i][1][1]))==2){
        fscanf(anglefile," ");
        fscanf(anglefile,"alpha");
        fscanf(anglefile," = ");
        fscanf(anglefile,"%f",angles3D[i][0]+1);
        fscanf(anglefile,"   ");
        fscanf(anglefile,"dim2=%1d%1d",&(combidim[i][0][2]),&(combidim[i][1][2]));
        fscanf(anglefile," ");
        fscanf(anglefile,"alpha");
        fscanf(anglefile," = ");
        fscanf(anglefile,"%f",angles3D[i][0]+2);
        fscanf(anglefile,"\n");
        i++;
      }
    }
    else {
      STOPMSG("not implemented code 8");
    }
  }
  totalnumproj=i;
  if(totalnumproj==0){
    STOPMSG("angles.dat is empty")
  }
  fclose(anglefile);
  if(fullrun){
    validanswer=0;
    i1=0;
    while(validanswer==0){
      GETINT("Choose an option:\n Measure until convergence (1)\n Measure a fixed number of projections (2)",i2)
      if(i2==1){
        validanswer=1;
        convergence=1;
      }
      else if(i2==2){
        validanswer=1;
	while((fixednumproj<=0)||(fixednumproj>totalnumproj)){
          GETINT("How many of these angle sets shall be used?",fixednumproj)
        }
      }
    }
  }
  else if(setuponly){
    while((fixednumproj<=0)||(fixednumproj>totalnumproj)){
      GETINT("How many of these angle sets shall be used?",fixednumproj)
    }
  }
}

/*---------------correct wrong angles-----------------*/

for(i=0;i<totalnumproj;i++){
  for(j=numofangles[1]-1;j>=0;j--){
    if(angles[i][j]<0.0) angles[i][j]=0.0;
    if(angles[i][j]>90.0) angles[i][j]=90.0;
    if((angles[i][j]==90.0) && (j>0)){
      for(k=j-1;k>=0;k--){
        angles[i][k]=0.0;
      }
      j=-1;
    }
  }
}

/*---------------analzsis of angles-----------------*/

for(i=0;i<totalnumproj;i++){
  ptomeasure[i]=1;
  nondangles[i]=0;
  if(dimproj==2){
    for(j=0;j<numofangles[1];j++){
      if((angles[i][j]>0.0) && (angles[i][j]<90.0)) {
        ptomeasure[i]*=2;
        nondangles[i]++;
        plusminus[i][j]=nondangles[i];
      }
      else plusminus[i][j]=0;
    }
    if(pow(2,nondangles[i])!=ptomeasure[i]){
      STOPMSG("bad programming 1")
    }
  }
  else if(dimproj==3){
    nondangles3D[i][1]=0;
    nondangles3D[i][2]=0;
    if(dimensionality==4){
      if((angles3D[i][0][1]>0.0) && (angles3D[i][0][1]<90.0)){
        ptomeasure[i]*=2;
        nondangles3D[i][1]++;
        plusminus3D[i][0][1]=nondangles3D[i][1];
      }
      else plusminus3D[i][0][1]=0;
    }
    else if(dimensionality==5){
      if((angles3D[i][0][1]>0.0) && (angles3D[i][0][1]<90.0)){
        ptomeasure[i]*=2;
        nondangles3D[i][1]++;
        plusminus3D[i][0][1]=nondangles3D[i][1];
      }
      else plusminus3D[i][0][1]=0;
      if((angles3D[i][0][2]>0.0) && (angles3D[i][0][2]<90.0)){
        ptomeasure[i]*=2;
        nondangles3D[i][2]++;
        plusminus3D[i][0][2]=nondangles3D[i][2];
      }
      else plusminus3D[i][0][2]=0;
    }
    else {
      STOPMSG("not implemented code 17");
    }
    if(pow(2,nondangles3D[i][1]+nondangles3D[i][2])!=ptomeasure[i]){
      STOPMSG("bad programming 2")
    }
  }
}

if(dimproj==3){
  if(dimensionality==4){
    for(i=0;i<totalnumproj;i++){
      for(j=0;j<dimensionality-1;j++){
    	dimused3D[i][j][1]=0;
    	dimused3D[i][j][2]=0;
      }
      if(angles3D[i][0][1]==90.0){
        dimused3D[i][combidim[i][0][1]-1][1]=1;
      }
      else if(angles3D[i][0][1]>0.0){
        dimused3D[i][combidim[i][0][1]-1][1]=1;
        dimused3D[i][combidim[i][1][1]-1][1]=1;
      }
      else{
        dimused3D[i][combidim[i][1][1]-1][1]=1;
      }
      k=0;
      for(j=0;j<dimensionality-1;j++){
        k+=dimused3D[i][j][1];
      }
      if(k!=nondangles3D[i][1]+1){
        sprintf(errorstring,"i=%d k=%d nondangles=%d angle=%f",i+1,k,nondangles3D[i][1],angles3D[i][0][1]);
        STOPMSG(errorstring)
      }
      dimused3D[i][combidim[i][0][2]-1][2]=1;
      k=0;
      for(j=0;j<dimensionality-1;j++){
        k+=dimused3D[i][j][2];
      }
      if(k!=nondangles3D[i][2]+1){
        sprintf(errorstring,"i=%d k=%d nondangles=%d angle=%f",i+1,k,nondangles3D[i][2],angles3D[i][0][2]);
        STOPMSG(errorstring)
      }
    }
  }
  else if(dimensionality==5){
    for(i=0;i<totalnumproj;i++){
      for(j=0;j<dimensionality-1;j++){
    	dimused3D[i][j][1]=0;
    	dimused3D[i][j][2]=0;
      }
      for(kk=1;kk<=2;kk++){
    	if(angles3D[i][0][kk]==90.0){
    	  dimused3D[i][combidim[i][0][kk]-1][kk]=1;
    	}
    	else if(angles3D[i][0][kk]>0.0){
    	  dimused3D[i][combidim[i][0][kk]-1][kk]=1;
    	  dimused3D[i][combidim[i][1][kk]-1][kk]=1;
    	}
    	else{
    	  dimused3D[i][combidim[i][1][kk]-1][kk]=1;
    	}
    	k=0;
    	for(j=0;j<dimensionality-1;j++){
    	  k+=dimused3D[i][j][kk];
    	}
    	if(k!=nondangles3D[i][kk]+1){
    	  sprintf(errorstring,"i=%d kk=%d k=%d nondangles=%d angle=%f",i+1,kk,k,nondangles3D[i][kk],angles3D[i][0][kk]);
    	  STOPMSG(errorstring)
    	}
      }
    }
  }
  else{
    STOPMSG("not implemented code 18");
  }
  for(i=0;i<totalnumproj;i++){
    for(j=0;j<MAXDIMENSIONS;j++){
      dimusedlett3D[i][1][j]='\0';
      dimusedlett3D[i][2][j]='\0';
    }
    for(kk=1;kk<=2;kk++){
      j=0;
      k=0;
      while(j<nondangles3D[i][kk]+1){
  	if(dimused3D[i][k][kk]){
  	  dimusedlett3D[i][kk][j]='y'-dimensionality+2+k;
  	  j++;
  	}
  	k++;
  	if(k>1000){
  	  STOPMSG("program error");
  	}
      }
    }
  }
}
else{
  for(i=0;i<totalnumproj;i++){
    for(j=0;j<dimensionality-1;j++){
      dimused[i][j]=0;
    }
    for(j=numofangles[1]-1;j>=0;j--){
      if(angles[i][j]==90.0){
  	dimused[i][dimensionality-3-j]=1;
  	j=-1;
      }
      else if(angles[i][j]>0.0){
  	dimused[i][dimensionality-3-j]=1;
      }
      if((j==0) && (angles[i][j]<90.0)){
  	dimused[i][dimensionality-2]=1;
      }
    }
    k=0;
    for(j=0;j<dimensionality-1;j++){
      k+=dimused[i][j];
    }
    if(k!=nondangles[i]+1){
      STOPMSG("bad programming 4")
    }
  }
  for(i=0;i<totalnumproj;i++){
    for(j=0;j<MAXDIMENSIONS;j++){
      dimusedlett[i][j]='\0';
    }
    j=0;
    k=0;
    while(j<nondangles[i]+1){
      if(dimused[i][k]){
    	dimusedlett[i][j]='y'-dimensionality+2+k;
    	j++;
      }
      k++;
      if(k>1000){
    	STOPMSG("program error");
      }
    }
  }
}

if((logfile=fopen(log,"a"))==NULL){
  STOPMSG("can not open logfile\n");
}
fprintf(logfile,"--------------------------------------------------\n");
fprintf(logfile,"!               APSY log file                    !\n");
fprintf(logfile,"--------------------------------------------------\n");
if((convergence) && (totalnumproj==1)) fprintf(logfile,"The following set of angles is used:\n\n");
else if (convergence) fprintf(logfile,"\nThe following %d sets of angles are used:\n\n",totalnumproj);
else if (fixednumproj==1) fprintf(logfile,"The following set of angles is used:\n\n");
else fprintf(logfile,"\nThe following %d sets of angles are used:\n\n",fixednumproj);
if(convergence) limit=totalnumproj;
else limit=fixednumproj;
dataset=expno;
for(j=0;j<limit;j++){
  fprintf(logfile,"Round %2d:",j+1);
  if(dimproj==2){
    for(i=0;i<numofangles[1];i++) {
      fprintf(logfile," %s = %4.1f  ",greekanglenames[i+1],angles[j][i]);
    }
  }
  else{
    if(dimensionality==4){
      fprintf(logfile," dim1=%d%d",combidim[j][0][1],combidim[j][1][1]);
      fprintf(logfile," alpha = %4.1f  ",angles3D[j][0][1]);
      fprintf(logfile," dim2=%d",combidim[j][0][2]);
    }
    else if(dimensionality==5){
      fprintf(logfile," dim1=%d%d",combidim[j][0][1],combidim[j][1][1]);
      fprintf(logfile," alpha = %4.1f  ",angles3D[j][0][1]);
      fprintf(logfile," dim2=%d%d",combidim[j][0][2],combidim[j][1][2]);
      fprintf(logfile," alpha = %4.1f  ",angles3D[j][0][2]);
    }
    else {
      STOPMSG("not implemented code 22");
    }
  }
  if(ptomeasure[j]==1) fprintf(logfile,"  1 spectrum in expno\t%d.\t", dataset);
  else fprintf(logfile,"  %d spectra in expnos\t%d-%d.\t",ptomeasure[j], dataset, dataset+ptomeasure[j]-1);
  fprintf(logfile," Dims: ");
  if(dimproj==2){
    for(i=0;i<dimensionality-1;i++) {
      fprintf(logfile,"%d",dimused[j][i]);
    }
  }
  else {
    for(i=0;i<dimensionality-1;i++) {
      fprintf(logfile,"%d",dimused3D[j][i][1]);
    }
    fprintf(logfile," and ");
    for(i=0;i<dimensionality-1;i++) {
      fprintf(logfile,"%d",dimused3D[j][i][2]);
    }
  }
  if(dimproj==2){
    fprintf(logfile," %s %d",dimusedlett[j],nondangles[j]);
  }
  else {
    fprintf(logfile," %s %d",dimusedlett3D[j][1],nondangles3D[j][1]);
    fprintf(logfile," %s %d",dimusedlett3D[j][2],nondangles3D[j][2]);
  }
  fprintf(logfile,"\n");
  dataset+=ptomeasure[j];
}
fprintf(logfile,"\n\n");
fclose(logfile);

while(apsycomplete==0){
  transient++;
/*------------------------------------generate the next dataset----------------*/
  if(transient>totalnumproj){
    STOPMSG("Not enough angles given for convergence.\n")
  }
  if(transient>1) {
    IEXPNO
  }
  basexpno=expno;
  VIEWDATA
  STOREPAR1("FnMODE",5)
  trosy=0;
  if(dimproj==2){
    if(strncmp(PLATFORM,"TOPSPIN",7)==0){
      if(dimensionality==3){
        FETCHPAR("CNST 41",&sweepx)
        FETCHPAR("CNST 42",&sweepy)
      }
      else if(dimensionality==4){
        FETCHPAR("CNST 41",&sweepw)
        FETCHPAR("CNST 42",&sweepx)
        FETCHPAR("CNST 43",&sweepy)
      }
      else if(dimensionality==5){
        FETCHPAR("CNST 41",&sweepv)
        FETCHPAR("CNST 42",&sweepw)
        FETCHPAR("CNST 43",&sweepx)
        FETCHPAR("CNST 44",&sweepy)
      }
      else if(dimensionality==6){
        FETCHPAR("CNST 41",&sweepu)
        FETCHPAR("CNST 42",&sweepv)
        FETCHPAR("CNST 43",&sweepw)
        FETCHPAR("CNST 44",&sweepx)
        FETCHPAR("CNST 45",&sweepy)
      }
      else if(dimensionality==7){
        FETCHPAR("CNST 41",&sweept)
        FETCHPAR("CNST 42",&sweepu)
        FETCHPAR("CNST 43",&sweepv)
        FETCHPAR("CNST 44",&sweepw)
        FETCHPAR("CNST 45",&sweepx)
        FETCHPAR("CNST 46",&sweepy)
      }
    }
    else if(strcmp(PLATFORM,"XWINNMR 3.5")==0){
      if(dimensionality==3){
  	FETCHPAR("CNST 21",&sweepx)
  	FETCHPAR("CNST 22",&sweepy)
      }
      else if(dimensionality==4){
  	FETCHPAR("CNST 21",&sweepw)
  	FETCHPAR("CNST 22",&sweepx)
  	FETCHPAR("CNST 23",&sweepy)
      }
      else if(dimensionality==5){
  	FETCHPAR("CNST 21",&sweepv)
  	FETCHPAR("CNST 22",&sweepw)
  	FETCHPAR("CNST 23",&sweepx)
  	FETCHPAR("CNST 24",&sweepy)
      }
      else if(dimensionality==6){
  	FETCHPAR("CNST 21",&sweepu)
  	FETCHPAR("CNST 22",&sweepv)
  	FETCHPAR("CNST 23",&sweepw)
  	FETCHPAR("CNST 24",&sweepx)
  	FETCHPAR("CNST 25",&sweepy)
      }
      else if(dimensionality==7){
  	FETCHPAR("CNST 21",&sweept)
  	FETCHPAR("CNST 22",&sweepu)
  	FETCHPAR("CNST 23",&sweepv)
  	FETCHPAR("CNST 24",&sweepw)
  	FETCHPAR("CNST 25",&sweepx)
  	FETCHPAR("CNST 26",&sweepy)
      }
    }
    else {
      STOPMSG("Platform not defined")
    }
    if(dimensionality==3){
      sweep1 =pow(sin(0.017453*angles[transient-1][0])*sweepx,2);
      sweep1+=pow(cos(0.017453*angles[transient-1][0])*sweepy,2);
      STOREPAR("CNST 1",angles[transient-1][0])
    }
    else if(dimensionality==4){
      sweep1 =pow(sin(0.017453*angles[transient-1][1])*sweepw,2);
      sweep1+=pow(sin(0.017453*angles[transient-1][0])*cos(0.017453*angles[transient-1][1])*sweepx,2);
      sweep1+=pow(cos(0.017453*angles[transient-1][0])*cos(0.017453*angles[transient-1][1])*sweepy,2);
      STOREPAR("CNST 1",angles[transient-1][0])
      STOREPAR("CNST 2",angles[transient-1][1])
    }
    else if(dimensionality==5){
      sweep1 =pow(sin(0.017453*angles[transient-1][2])*sweepv,2);
      sweep1+=pow(sin(0.017453*angles[transient-1][1])*cos(0.017453*angles[transient-1][2])*sweepw,2);
      sweep1+=pow(sin(0.017453*angles[transient-1][0])*cos(0.017453*angles[transient-1][1])*cos(0.017453*angles[transient-1][2])*sweepx,2);
      sweep1+=pow(cos(0.017453*angles[transient-1][0])*cos(0.017453*angles[transient-1][1])*cos(0.017453*angles[transient-1][2])*sweepy,2);
      STOREPAR("CNST 1",angles[transient-1][0])
      STOREPAR("CNST 2",angles[transient-1][1])
      STOREPAR("CNST 3",angles[transient-1][2])
    }
    else if(dimensionality==6){
      sweep1 =pow(sin(0.017453*angles[transient-1][3])*sweepu,2);
      sweep1+=pow(sin(0.017453*angles[transient-1][2])*cos(0.017453*angles[transient-1][3])*sweepv,2);
      sweep1+=pow(sin(0.017453*angles[transient-1][1])*cos(0.017453*angles[transient-1][2])*cos(0.017453*angles[transient-1][3])*sweepw,2);
      sweep1+=pow(sin(0.017453*angles[transient-1][0])*cos(0.017453*angles[transient-1][1])*cos(0.017453*angles[transient-1][2])*cos(0.017453*angles[transient-1][3])*sweepx,2);
      sweep1+=pow(cos(0.017453*angles[transient-1][0])*cos(0.017453*angles[transient-1][1])*cos(0.017453*angles[transient-1][2])*cos(0.017453*angles[transient-1][3])*sweepy,2);    
      STOREPAR("CNST 1",angles[transient-1][0])
      STOREPAR("CNST 2",angles[transient-1][1])
      STOREPAR("CNST 3",angles[transient-1][2])
      STOREPAR("CNST 4",angles[transient-1][3])
    }
    else if(dimensionality==7){
      sweep1 =pow(sin(0.017453*angles[transient-1][4])*sweept,2);
      sweep1+=pow(sin(0.017453*angles[transient-1][3])*cos(0.017453*angles[transient-1][4])*sweepu,2);
      sweep1+=pow(sin(0.017453*angles[transient-1][2])*cos(0.017453*angles[transient-1][3])*cos(0.017453*angles[transient-1][4])*sweepv,2);
      sweep1+=pow(sin(0.017453*angles[transient-1][1])*cos(0.017453*angles[transient-1][2])*cos(0.017453*angles[transient-1][3])*cos(0.017453*angles[transient-1][4])*sweepw,2);
      sweep1+=pow(sin(0.017453*angles[transient-1][0])*cos(0.017453*angles[transient-1][1])*cos(0.017453*angles[transient-1][2])*cos(0.017453*angles[transient-1][3])*cos(0.017453*angles[transient-1][4])*sweepx,2);
      sweep1+=pow(cos(0.017453*angles[transient-1][0])*cos(0.017453*angles[transient-1][1])*cos(0.017453*angles[transient-1][2])*cos(0.017453*angles[transient-1][3])*cos(0.017453*angles[transient-1][4])*sweepy,2);    
      STOREPAR("CNST 1",angles[transient-1][0])
      STOREPAR("CNST 2",angles[transient-1][1])
      STOREPAR("CNST 3",angles[transient-1][2])
      STOREPAR("CNST 4",angles[transient-1][3])
      STOREPAR("CNST 5",angles[transient-1][4])
    }
  }
  else if(dimproj==3){
    if(strcmp(PLATFORM,"TOPSPIN",7)==0){
      if(dimensionality==4){
        FETCHPAR("CNST 41",&(sweep_a[1]))
        FETCHPAR("CNST 42",&(sweep_a[2]))
        FETCHPAR("CNST 43",&(sweep_a[3]))
      }
      else if(dimensionality==5){
        FETCHPAR("CNST 41",&(sweep_a[1]))
        FETCHPAR("CNST 42",&(sweep_a[2]))
        FETCHPAR("CNST 43",&(sweep_a[3]))
        FETCHPAR("CNST 44",&(sweep_a[4]))
      }
      else{
        STOPMSG("dimensionality not implemented code 26");
      }
    }
    else if(strcmp(PLATFORM,"XWINNMR 3.5")==0){
      if(dimensionality==4){
        FETCHPAR("CNST 21",&(sweep_a[1]))
        FETCHPAR("CNST 22",&(sweep_a[2]))
        FETCHPAR("CNST 23",&(sweep_a[3]))
      }
      else if(dimensionality==5){
        FETCHPAR("CNST 21",&(sweep_a[1]))
        FETCHPAR("CNST 22",&(sweep_a[2]))
        FETCHPAR("CNST 23",&(sweep_a[3]))
        FETCHPAR("CNST 24",&(sweep_a[4]))
      }
      else{
        STOPMSG("dimensionality not implemented code 26");
      }
    }
    else {
      STOPMSG("Platform not defined")
    }
    if(dimensionality==4){
      sweep1 =pow(sin(0.017453*angles3D[transient-1][0][1])*sweep_a[combidim[transient-1][0][1]],2);
      sweep1+=pow(cos(0.017453*angles3D[transient-1][0][1])*sweep_a[combidim[transient-1][1][1]],2);
      sweep3 =pow(sweep_a[combidim[transient-1][0][2]],2);
      STOREPAR("CNST 1", angles3D[transient-1][0][1])
    }
    else if(dimensionality==5){
      sweep1 =pow(sin(0.017453*angles3D[transient-1][0][1])*sweep_a[combidim[transient-1][0][1]],2);
      sweep1+=pow(cos(0.017453*angles3D[transient-1][0][1])*sweep_a[combidim[transient-1][1][1]],2);
      sweep3 =pow(sin(0.017453*angles3D[transient-1][0][2])*sweep_a[combidim[transient-1][0][2]],2);
      sweep3+=pow(cos(0.017453*angles3D[transient-1][0][2])*sweep_a[combidim[transient-1][1][2]],2);
      STOREPAR("CNST 1", angles3D[transient-1][0][1])
      STOREPAR("CNST 2", angles3D[transient-1][0][2])
    }
  }
  sweep1=sqrt(sweep1);
  STOREPAR1("SWH",sweep1)
  STOREPAR("IN 1",1.0/sweep1);
  STOREPAR1("OFFSET",0.0)
  if(dimproj==3) {
    sweep3=sqrt(sweep3);
    STOREPAR3("SWH",sweep3)
    STOREPAR3("OFFSET",0.0)
  }
  if((trosydim!=0) && (dimused[transient-1][dimensionality-2]==1)){
    trosy=1;
    STOREPAR1("FnMODE",6)
  }
  numofp=ptomeasure[transient-1];
/*----------------------------------creation of ZGoptns-------------------------------------------*/
  for(k=1;k<=numofp;k++){
    if(k>1){
      IEXPNO
    }
    if(dimproj==2){
      for(i=0;i<=30;i++){
        phasesused[i]='0';
      }
      sprintf(zgoptns,"-DP%s",dimusedlett[transient-1]+nondangles[transient-1]);
      for(j=0;j<=nondangles[transient-1];j++){
        if(tppi[k-1][5-nondangles[transient-1]+j]==1) phasesused[dimusedlett[transient-1][j]-'y'+dimensionality-2]='1';
        if(phaseoverlap[dimusedlett[transient-1][j]-'y'+dimensionality-2]=='0'){
          zgadd[3]=dimusedlett[transient-1][j]+'A'-'a';
          if(tppi[k-1][5-nondangles[transient-1]+j]==1) strcat(zgoptns,zgadd);
        }
      }
      if(poexists){
        j=0;
        i=0;
        strcat(zgoptns," -D");
        while(j<poexists){
          if(phaseoverlap[i]=='1'){
    	    j++;
    	    zgadd2[0]='Y'-dimensionality+2+i;
    	    strcat(zgoptns,zgadd2);
    	    zgadd2[0]=phasesused[i];
    	    strcat(zgoptns,zgadd2);
    	  }
          i++;
        }
      }
      if(receiver[k-1]==2) strcat(zgoptns," -DA2");
    }
    else{
      sprintf(zgoptns,"-D%s_%s",dimusedlett3D[transient-1][1],dimusedlett3D[transient-1][2]);
      if(dimensionality==4){
        if     ((combidim[transient-1][0][1]==1)&&(combidim[transient-1][1][1]==2)&&(combidim[transient-1][0][2]==3)) strcat(zgoptns," -DA");
        else if((combidim[transient-1][0][1]==1)&&(combidim[transient-1][1][1]==3)&&(combidim[transient-1][0][2]==2)) strcat(zgoptns," -DB");
        else if((combidim[transient-1][0][1]==2)&&(combidim[transient-1][1][1]==3)&&(combidim[transient-1][0][2]==1)) strcat(zgoptns," -DC");
	else {
	  STOPMSG("wrong usage of dimensions");
	}
      }
      else if(dimensionality==5){
        if(combidim[transient-1][0][1]!=1){
	  STOPMSG("wrong usage of dimensions");
	}
        if     ((combidim[transient-1][1][1]==2)&&(combidim[transient-1][0][2]==3)&&(combidim[transient-1][1][2]==4)) strcat(zgoptns," -DA");
        else if((combidim[transient-1][1][1]==3)&&(combidim[transient-1][0][2]==2)&&(combidim[transient-1][1][2]==4)) strcat(zgoptns," -DB");
        else if((combidim[transient-1][1][1]==4)&&(combidim[transient-1][0][2]==2)&&(combidim[transient-1][1][2]==3)) strcat(zgoptns," -DC");
	else {
	  STOPMSG("wrong usage of dimensions");
	}
      }
      else {
        STOPMSG("not implemented code 34");
      }
      for(kk=1;kk<=2;kk++){
        for(j=0;j<=nondangles3D[transient-1][kk];j++){
          zgadd[3]=dimusedlett3D[transient-1][kk][j]+'A'-'a';
	  if(numofp==2){
            if((k==2)&&(nondangles3D[transient-1][kk]==1)) strcat(zgoptns,zgadd);
	  }
	  else if(numofp==4){
            if((k==2)&&(kk==1)) strcat(zgoptns,zgadd);
            if((k==3)&&(kk==2)) strcat(zgoptns,zgadd);
            if(k==4) strcat(zgoptns,zgadd);
	  }
	  else if(numofp!=1){
	    STOPMSG("not possible");
	  }
	}
      }
      if((k==2)||(k==3)) strcat(zgoptns," -DA2");
    }
    STOREPAR("ZGOPTNS",zgoptns)
  }
/*-----------------------------------write PROSA file----------------------------------------------*/

  sprintf(prosa,"%s/data/%s/nmr/%s/%d/apsy.pro",disk,user,name,basexpno);
  if((prosafile=fopen(prosa, "w"))==NULL) {
    STOPMSG("can't create file apsy.pro to write\n");
  }
  if(trosy==1){
    for(countter=1;countter<=numofp;countter++){
      fprintf(prosafile,"read %s ../%d/ser %dc %dc\n",byteorder,basexpno-1+countter,pts2complex,pts1complex);
      fprintf(prosafile,"dimension 2\n");
      fprintf(prosafile,"phase 90\n");
      fprintf(prosafile,"dimension 1\n");
      fprintf(prosafile,"read %s ../%d/ser combine 1.0 1.0\n",byteorder,basexpno-1+countter);
      fprintf(prosafile,"select * i\n");
      fprintf(prosafile,"phase 90\n");
      fprintf(prosafile,"select\n");
      fprintf(prosafile,"phase -90\n");
      fprintf(prosafile,"write %s ../%d/Tser%d\n\n",byteorder,basexpno,countter);
    }
    for(countter=1;countter<=numofp;countter++){
      fprintf(prosafile,"read %s ../%d/Tser1 %dc %dc\n",byteorder,basexpno,pts2complex,pts1complex);
      for(count2=2;count2<=numofp;count2++){
	fprintf(prosafile,"read %s ../%d/Tser%d combine 1 %d\n",byteorder,basexpno,count2,addtheorem[countter-1][count2-1]);
      }
      fprintf(prosafile,"write %s ../%d/ser%d\n\n",byteorder,basexpno,countter);
    }
  }
  else{
    for(countter=1;countter<=numofp;countter++){
      fprintf(prosafile,"read %s ../%d/ser %dc %dc",byteorder,basexpno,pts2complex,pts1complex);
      if(dimproj==3) fprintf(prosafile," %dc",pts3complex);
      fprintf(prosafile,"\n");
      for(count2=2;count2<=numofp;count2++){
	fprintf(prosafile,"read %s ../%d/ser combine 1 %d\n",byteorder,basexpno-1+count2,addtheorem[countter-1][count2-1]);
      }
      fprintf(prosafile,"write %s ../%d/ser%d\n\n",byteorder,basexpno,countter);
    }
  }
  for(countter=1;countter<=numofp;countter++){
    fprintf(prosafile,"read %s ../%d/ser%d %dc %dc",byteorder,basexpno,countter,pts2complex,pts1complex);
    if(dimproj==3) fprintf(prosafile," %dc",pts3complex);
    fprintf(prosafile,"\n");
    fprintf(prosafile,"decim      = %d\n",decimm);
    fprintf(prosafile,"gdp        = %6.3f\n",groupdelay);
    fprintf(prosafile,"#DIMENSION: 1\n");
    fprintf(prosafile,"si(1)      = %d\n",si2);
    fprintf(prosafile,"aqm(1)    := RSH\n");
    fprintf(prosafile,"w0(1)	= 1\n");
    fprintf(prosafile,"delta(1)   = %12.10f\n", 1/SWhertz);
    fprintf(prosafile,"ppmmax(1)  = %8.3f\n", 0.5*SWhertz);
    fprintf(prosafile,"#DIMENSION: 2\n");
    fprintf(prosafile,"si(2)      = $ptindi\n");
    fprintf(prosafile,"aqm(2)    := ST/TP\n");
    fprintf(prosafile,"w0(2)	  = 1\n");
    fprintf(prosafile,"delta(2)   = %12.10f\n", 1/sweep1);
    fprintf(prosafile,"ppmmax(2)  =  (0.5/$delta(2))*(1.0+1.0/($si(2)*1.0))\n");
    if(dimproj==3){
      fprintf(prosafile,"#DIMENSION: 3\n");
      fprintf(prosafile,"si(3)      = $ptindi3\n");
      fprintf(prosafile,"aqm(3)    := ST/TP\n");
      fprintf(prosafile,"w0(3)	  = 1\n");
      fprintf(prosafile,"delta(3)   = %12.10f\n", 1/sweep3);
      fprintf(prosafile,"ppmmax(3)  =  (0.5/$delta(3))*(1.0+1.0/($si(3)*1.0))\n");
    }
    fprintf(prosafile,"ns	      = %d\n",numberscans);
    fprintf(prosafile,"set dim_order= 1");
    if(dimproj==3) fprintf(prosafile," 2 3");
    fprintf(prosafile,"\n");
    fprintf(prosafile,"set n_dim= %d\n",dimproj);
    fprintf(prosafile,"dimension 1\n");
    fprintf(prosafile,"multiply 0.5 1\n");
    fprintf(prosafile,"window sin 75 %6.3f\n",groupdelay);
    fprintf(prosafile,"ft %d %d %d $phcor %6.1f full\n", si2, 1+stsr2, stsr2+finalsi2, pphhase2);
    fprintf(prosafile,"dimension 2\n");
    if(trosy==0){
      fprintf(prosafile,"multiply -1 2 %d 2\n",pts1complex);
    }
    fprintf(prosafile,"if ('$linpred'=='y') predict lpsvd $shapes $ptspredict\n");
    fprintf(prosafile,"if ('$mirrlp'=='y') mirrorlp $shapes $ptspredict 2\n");
    fprintf(prosafile,"multiply 0.5 1\n");
    fprintf(prosafile,"window $wfunction\n");
    fprintf(prosafile,"ft $si(2) 1 $si(2) 0 0 full\n");
    if(dimproj==3){
      fprintf(prosafile,"dimension 3\n");
      fprintf(prosafile,"multiply -1 2 %d 2\n",pts3complex);
      fprintf(prosafile,"multiply 0.5 1\n");
      fprintf(prosafile,"window sin 75\n");
      fprintf(prosafile,"ft $si(3) 1 $si(3) 0 0 full\n");
      fprintf(prosafile,"dimension 3\n");
      fprintf(prosafile,"if ('$baseline3'=='y') cflatt flatt 10 5 polynom 2\n");
      fprintf(prosafile,"dimension 2\n");
    }
    fprintf(prosafile,"if ('$baseline2'=='y') cflatt flatt 10 5 polynom 2\n");
    fprintf(prosafile,"dimension 1\n");
    fprintf(prosafile,"cflatt iterative cftw 5 $si(1)\n");
    if(dimproj==3) fprintf(prosafile,"dimension 1 3 2\n");
    fprintf(prosafile,"scalf = 100/snoise\n");
    fprintf(prosafile,"multiply $scalf\n");
    fprintf(prosafile,"write easy16 /nmr/data/%s/%s/$folder/spec%d byteswap\n\n",user,name,basexpno+countter-1);
  }
  fclose(prosafile);
  if((megaprosa=fopen(megaprosname,"a"))==NULL) {
    STOPMSG("can not open megaprosa\n")
  }
  fprintf(megaprosa,"../%d/apsy\n",basexpno);
  fclose(megaprosa);
  
  VIEWDATA
/*--------------------------run experiments----------------------------------*/
  if(fullrun){
     if((logfile=fopen(log,"a"))==NULL){
       STOPMSG("can not open logfile\n");
     }
     fprintf(logfile,"\nNMR experiments started");
     for(i=0;i<numofangles[1];i++) {
       fprintf(logfile,",  %s = ",greekanglenames[i+1]);
       if(plusminus[transient-1][i]!=0) fprintf(logfile,"+/-");
       fprintf(logfile,"%4.1f",angles[transient-1][i]);
     }
     fprintf(logfile,"\n");
     fclose(logfile);
     for(i=1;i<numofp;i++){
       DEXPNO
     }
     if((logfile=fopen(log,"a"))==NULL){
       STOPMSG("can not open logfile\n");
     }
     fprintf(logfile,"Experiment 1 out of %d started.\n",numofp);
     fclose(logfile);
     ZG
     for(i=1;i<numofp;i++){
       IEXPNO
       if((logfile=fopen(log,"a"))==NULL){
         STOPMSG("can not open logfile\n");
       }
       fprintf(logfile,"Experiment %d out of %d started.\n",i+1,numofp);
       fclose(logfile);
       ZG
     }
     if((logfile=fopen(log,"a"))==NULL){
       STOPMSG("can not open logfile\n");
     }
     fprintf(logfile,"Measurements ended.\n");
     fclose(logfile);
  }
  for(i=1;i<numofp;i++){
    DEXPNO
  }
  if(fullrun){
  /*process the just recorded or the current experiments*/
    if((logfile=fopen(log,"a"))==NULL){
      STOPMSG("can not open logfile\n");
    }
    fprintf(logfile,"\nSpectra processing started.\n");
    fclose(logfile);
    for(j=1;j<numofp;j++){
      XF2
      IEXPNO
    }
    XF2
    for(j=1;j<numofp;j++){
      DEXPNO
    }
    if(numofp==2){
      STOREPAR("ALPHA",1.0)
      STOREPAR("GAMMA",1.0)
      DATASET2(name, basexpno+1, 1, disk, user)
      ADD2D
      IEXPNO
      STOREPAR("ALPHA",-2.0)
      STOREPAR("GAMMA",1.0)
      DATASET2(name, basexpno, 1, disk, user)
      ADD2D
    }
    else if(numofp==4){
      STOREPAR("ALPHA",1.0)
      STOREPAR("GAMMA",1.0)
      DATASET2(name, basexpno+1, 1, disk, user)
      ADD2D
      IEXPNO
      STOREPAR("ALPHA",-2.0)
      STOREPAR("GAMMA",1.0)
      DATASET2(name, basexpno, 1, disk, user)
      ADD2D
      IEXPNO
      STOREPAR("ALPHA",1.0)
      STOREPAR("GAMMA",1.0)
      DATASET2(name, basexpno+3, 1, disk, user)
      ADD2D
      IEXPNO
      STOREPAR("ALPHA",-2.0)
      STOREPAR("GAMMA",1.0)
      DATASET2(name, basexpno+2, 1, disk, user)
      ADD2D
      DEXPNO
      DEXPNO
      DEXPNO
      STOREPAR("ALPHA",1.0)
      STOREPAR("GAMMA",1.0)
      DATASET2(name, basexpno+2, 1, disk, user)
      ADD2D
      IEXPNO
      STOREPAR("ALPHA",1.0)
      STOREPAR("GAMMA",1.0)
      DATASET2(name, basexpno+3, 1, disk, user)
      ADD2D
      IEXPNO
      STOREPAR("ALPHA",-2.0)
      STOREPAR("GAMMA",1.0)
      DATASET2(name, basexpno, 1, disk, user)
      ADD2D
      IEXPNO
      STOREPAR("ALPHA",-2.0)
      STOREPAR("GAMMA",1.0)
      DATASET2(name, basexpno+1, 1, disk, user)
      ADD2D
    }
    else if(numofp==8){
      STOREPAR("ALPHA",1.0)
      STOREPAR("GAMMA",1.0)
      DATASET2(name, basexpno+1, 1, disk, user)
      ADD2D
      IEXPNO
      STOREPAR("ALPHA",1.0)
      STOREPAR("GAMMA",1.0)
      DATASET2(name, basexpno+1, 1, disk, user)
      ADD2D
      IEXPNO
      STOREPAR("ALPHA",1.0)
      STOREPAR("GAMMA",1.0)
      DATASET2(name, basexpno+1, 1, disk, user)
      ADD2D
      IEXPNO
      STOREPAR("ALPHA",1.0)
      STOREPAR("GAMMA",1.0)
      DATASET2(name, basexpno+1, 1, disk, user)
      ADD2D
      IEXPNO
      STOREPAR("ALPHA",1.0)
      STOREPAR("GAMMA",1.0)
      DATASET2(name, basexpno+1, 1, disk, user)
      ADD2D
      IEXPNO
      STOREPAR("ALPHA",-2.0)
      STOREPAR("GAMMA",1.0)
      DATASET2(name, basexpno, 1, disk, user)
      ADD2D
      IEXPNO
      STOREPAR("ALPHA",1.0)
      STOREPAR("GAMMA",1.0)
      DATASET2(name, basexpno+3, 1, disk, user)
      ADD2D
      IEXPNO
      STOREPAR("ALPHA",-2.0)
      STOREPAR("GAMMA",1.0)
      DATASET2(name, basexpno+2, 1, disk, user)
      ADD2D
      DEXPNO
      DEXPNO
      DEXPNO
      DEXPNO
      DEXPNO
      DEXPNO
      DEXPNO
      STOREPAR("ALPHA",1.0)
      STOREPAR("GAMMA",1.0)
      DATASET2(name, basexpno+2, 1, disk, user)
      ADD2D
      IEXPNO
      STOREPAR("ALPHA",1.0)
      STOREPAR("GAMMA",1.0)
      DATASET2(name, basexpno+2, 1, disk, user)
      ADD2D
      IEXPNO
      STOREPAR("ALPHA",1.0)
      STOREPAR("GAMMA",1.0)
      DATASET2(name, basexpno+2, 1, disk, user)
      ADD2D
      IEXPNO
      STOREPAR("ALPHA",1.0)
      STOREPAR("GAMMA",1.0)
      DATASET2(name, basexpno+2, 1, disk, user)
      ADD2D
      IEXPNO
      STOREPAR("ALPHA",1.0)
      STOREPAR("GAMMA",1.0)
      DATASET2(name, basexpno+2, 1, disk, user)
      ADD2D
      IEXPNO
      STOREPAR("ALPHA",1.0)
      STOREPAR("GAMMA",1.0)
      DATASET2(name, basexpno+3, 1, disk, user)
      ADD2D
      IEXPNO
      STOREPAR("ALPHA",-2.0)
      STOREPAR("GAMMA",1.0)
      DATASET2(name, basexpno, 1, disk, user)
      ADD2D
      IEXPNO
      STOREPAR("ALPHA",-2.0)
      STOREPAR("GAMMA",1.0)
      DATASET2(name, basexpno+1, 1, disk, user)
      ADD2D
    }
    else if(numofp!=1){
      STOPMSG("Your data is not supported\n")
    }
    XF1
    ABS2
    ABS1
    for(j=1;j<numofp;j++){
      DEXPNO
      XF1
      ABS2
      ABS1
    }
    VIEWDATA
    if((logfile=fopen(log,"a"))==NULL){
      STOPMSG("can not open logfile\n");
    }
    fprintf(logfile,"Spectra processing finished.\n\n");
    fclose(logfile);
  }
  if((spectrafile=fopen(spectra,"a"))==NULL){
    STOPMSG("can not open spectrafile\n");
  }
  fprintf(spectrafile,"# %d\n",numofp);
  fclose(spectrafile);
  if((spectrafile2=fopen(spectra2,"a"))==NULL){
    STOPMSG("can not open spectrafile 2\n");
  }
  fprintf(spectrafile2,"# %d\n",numofp);
  fclose(spectrafile2);
  for(i=0;i<numofp;i++){
    if(i!=0) {
      IEXPNO
    }
    sprintf(spectrum,"%s/data/%s/nmr/%s/%d/pdata/1/2rr",disk,user,name,expno);
    sprintf(spectrum2,"spec%d",expno);
    if(fullrun){
      if((logfile=fopen(log,"a"))==NULL){
        STOPMSG("can not open logfile\n");
      }
      fprintf(logfile,"\nPeak picking spectrum with angles");
      for(k=1;k<=numofangles[1];k++){
        fprintf(logfile," %s = ",greekanglenames[k]);
        if(minuss[i][plusminus[transient-1][k-1]]!=0) fprintf(logfile,"-");
        fprintf(logfile,"%4.1f",angles[transient-1][k-1]);
      }
      fprintf(logfile,"\n");
      fclose(logfile);
      sprintf(execute,"%s -PeakPicking %s %s",gapro,mainpath,spectrum);
      system(execute);
      finished=0;
      while(finished==0){
        ssleep(0.1);
        if((communfile=fopen(communication,"r"))!=NULL){
          finished=1;
	  fclose(communfile);
        }
        if((errorfile=fopen(errorfname,"r"))!=NULL){
          fclose(errorfile);
          sprintf(execute,"rm %s",errorfname);
          system(execute);
          ssleep(0.1);
          STOPMSG("Unusual termination of GAPRO");
        }
      }
      sprintf(execute,"rm %s",communication);
      system(execute);
      ssleep(0.1);
      numofprojsofar++;
    }
    if((spectrafile=fopen(spectra,"a"))==NULL){
      STOPMSG("can not open spectrafile\n");
    }
    if((spectrafile2=fopen(spectra2,"a"))==NULL){
      STOPMSG("can not open spectrafile 2\n");
    }
    fprintf(spectrafile,"%s",spectrum);
    fprintf(spectrafile2,"%s",spectrum2);
    if(dimproj==2){
      for(k=1;k<=numofangles[1];k++){
        fprintf(spectrafile," %s=",greekanglenames[k]);
        fprintf(spectrafile2," %s=",greekanglenames[k]);
        if(minuss[i][plusminus[transient-1][k-1]]!=0) {
	  fprintf(spectrafile,"-");
          fprintf(spectrafile2,"-");
	}
        fprintf(spectrafile,"%f",angles[transient-1][k-1]);
        fprintf(spectrafile2,"%f",angles[transient-1][k-1]);
      }
      fprintf(spectrafile,"\n");
      fprintf(spectrafile2,"\n");
    }
    else if(dimproj==3){
      if(dimensionality==4){
        for(k=1;k<=numofangles[1];k++){
          fprintf(spectrafile," %s=",greekanglenames[k]);
          fprintf(spectrafile2," %s=",greekanglenames[k]);
          if(minuss[i][plusminus3D[transient-1][k-1][1]]!=0) {
	    fprintf(spectrafile,"-");
            fprintf(spectrafile2,"-");
	  }
          fprintf(spectrafile,"%f",angles3D[transient-1][k-1][1]);
          fprintf(spectrafile2,"%f",angles3D[transient-1][k-1][1]);
        }
        fprintf(spectrafile,"\n");
        fprintf(spectrafile2,"\n");
      }
      else {
        STOPMSG("not implemented code 46");
      }
    }
    fclose(spectrafile);
    fclose(spectrafile2);
  }
  if(fullrun){
    sprintf(execute,"%s -GeometricAnalysis %s",gapro, mainpath);
    system(execute);
    finished=0;
    while(finished==0){
      ssleep(0.1);
      if((communfile=fopen(communication,"r"))!=NULL){
    	finished=1;
	fscanf(communfile,"%d",peaksfound+transient);
        fclose(communfile);
      }
      if((errorfile=fopen(errorfname,"r"))!=NULL){
        fclose(errorfile);
        sprintf(execute,"rm %s",errorfname);
        system(execute);
        ssleep(0.1);
	STOPMSG("Unusual termination of GAPRO");
      }
    }
    sprintf(execute,"rm %s",communication);
    system(execute);
    ssleep(0.1);
    if((logfile=fopen(log,"a"))==NULL){
      STOPMSG("can not open logfile\n");
    }
    fclose(logfile);
  }
  if(convergence){
    if((logfile=fopen(log,"a"))==NULL){
      STOPMSG("can not open logfile\n");
    }
    fprintf(logfile,"\n--------------------------------------------------------\n");
    if((peaksfound[transient]>0) && (peaksfound[transient]<=peaksfound[transient-1])) {
      apsycomplete=1;
      fprintf(logfile,"\nConvergence was reached in round %d. Calculating final peak lists.\n",transient);
    }
    else {
      fprintf(logfile,"\nConvergence not reached in round %d. Measuring next angle set.\n",transient);
    }
    fprintf(logfile,"\n--------------------------------------------------------\n");
    fclose(logfile);
  }
  if(fixednumproj){
    if(transient==fixednumproj) apsycomplete=1;
  }
  if(fullrun){
    if((logfile=fopen(log,"a"))==NULL){
      STOPMSG("can not open logfile\n");
    }
    fprintf(logfile,"\n-------------------------------\n");
    fprintf(logfile,"\nSummary:\n");
    fprintf(logfile,"\n-------------------------------\n\n");
    for(i=1;i<=transient;i++){
      fprintf(logfile,"Round %d: %d Peaks\n",i,peaksfound[i]);
    }
    if ((convergence)&&(apsycomplete)) fprintf(logfile,"Converged with round %d\n",transient);
    if (fixednumproj>1) fprintf(logfile,"Totally %d rounds were run\n",transient);
    if (fixednumproj==1) fprintf(logfile,"1 round was run\n");
    if ((convergence)&&(apsycomplete)) {
      fprintf(logfile,"\nThe final %d-dimensional peak list was written:\n",dimensionality);
      fprintf(logfile,"%s/result.peaks\n",mainpath);
    }
    fclose(logfile);
  }
  if((communfile=fopen(communication,"r"))!=NULL){
    fclose(communfile);
    STOPMSG("resistant term 2");
  }
}

if((logfile=fopen(log,"a"))==NULL){
  STOPMSG("can not open logfile\n");
}
fprintf(logfile,"\n-------------------------------\n");
fprintf(logfile,"\nAPSY run ended successfully\n");
fprintf(logfile,"\n-------------------------------\n");
fclose(logfile);
QUIT