/*
 Authors 
 Martin Schlather, schlather@math.uni-mannheim.de

 Definition of auxiliary correlation functions 

Note:
 * Never use the below functions directly, but only by the functions indicated 
   in RFsimu.h, since there is gno error check (e.g. initialization of RANDOM)
 * VARIANCE, SCALE are not used here 
 * definitions for the random coin method can be found in MPPFcts.cc
 * definitions for genuinely anisotropic or nonsta     tionary models are in
   SophisticatedModel.cc; hyper models also in Hypermodel.cc

 Copyright (C) 2001 -- 2003 Martin Schlather
 Copyright (C) 2004 -- 2004 Yindeng Jiang & Martin Schlather
 Copyright (C) 2005 -- 2014 Martin Schlather

This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License
as published by the Free Software Foundation; either version 3
of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.  
*/


#include <math.h>
#include <stdio.h> 
#include <R_ext/Lapack.h>

#include "RF.h"
#include "Covariance.h"
#include "variogramAndCo.h"
//#include <R_ext/Applic.h>
//#include <R_ext/Utils.h>     
//#include <R_ext/BLAS.h> 





// $
void kappaS(int i, cov_model *cov, int *nr, int *nc){
  switch(i) {
  case DVAR : case DSCALE :
    *nr = *nc = 1;
    break;
  case DANISO :
    *nr = cov->xdimown;
    *nc = SIZE_NOT_DETERMINED;
    break;
  case DAUSER :
    *nr = SIZE_NOT_DETERMINED;
    *nc = cov->xdimown;
    break;
  case DPROJ : 
    *nr = SIZE_NOT_DETERMINED;
    *nc = 1;
    break;
  default : *nr = *nc = -1; 
  }
}


// simple transformations as variance, scale, anisotropy matrix, etc.  
void Siso(double *x, cov_model *cov, double *v){


  //  printf("call Siso %s %s %s\n", NAME(cov), NAME(cov->calling),
  //	 isDollar(cov) ? NAME(cov->sub[0]) : "-");
  //PMI(cov);

  cov_model *next = cov->sub[DOLLAR_SUB];
  int i,
    vdim = cov->vdim[0],
    vdimSq = vdim * vdim;
  double y,
    *aniso=P(DANISO),
    *scale =P(DSCALE),
    var = P0(DVAR);
  
  //  printf("\nSiso\n");

  y = *x;
  if (aniso != NULL) y = fabs(y * aniso[0]);

  if (scale != NULL) 
    y = scale[0]>0.0 ? y / scale[0] 
      : (y == 0.0 && scale[0]==0.0) ? 0.0 : RF_INF;
      
  // letzteres darf nur passieren wenn dim = 1!!
  COV(&y, next, v);

  //printf("end isoS\n");

  // PMI(cov);
  //  PMI(cov->calling->calling);
  // print("$ %f %f %f %d\n", *x, var, v[0], vdimSq);

  for (i=0; i<vdimSq; i++) v[i] *= var; 

  // print("$ %f %f %f %d\n", *x, var, v[0], vdimSq);
  //  error("von wo wird aufgerufen, dass nan in var?? und loc:x not given?");

}
  

// simple transformations as variance, scale, anisotropy matrix, etc.  
void logSiso(double *x, cov_model *cov, double *v, double *sign){
  cov_model *next = cov->sub[DOLLAR_SUB];
  int i,
    vdim = cov->vdim[0],
    vdimSq = vdim * vdim;
  double y, 
    *aniso=P(DANISO),
    *scale =P(DSCALE),
    logvar = log(P0(DVAR));
  
  y = *x;
  if (aniso != NULL) y = fabs(y * aniso[0]);

  if (scale != NULL) 
    y = scale[0]>0.0 ? y / scale[0] 
      : (y == 0.0 && scale[0]==0.0) ? 0.0 : RF_INF;
      
  LOGCOV(&y, next, v, sign);
  for (i=0; i<vdimSq; i++) v[i] += logvar; 
}
 
void Sstat(double *x, cov_model *cov, double *v){
  logSstat(x, cov, v, NULL);
}

void logSstat(double *x, cov_model *cov, double *v, double *sign){
  cov_model *next = cov->sub[DOLLAR_SUB],
    *Aniso = cov->kappasub[DAUSER];
  double 
    var = P0(DVAR),
    *scale =P(DSCALE), 
    *aniso=P(DANISO);
  int i,
    nproj = cov->nrow[DPROJ],
    vdim = cov->vdim[0],
    vdimSq = vdim * vdim;
 
  if (nproj > 0) {
    int *proj = PINT(DPROJ);
    ALLOC_DOLLAR(z, nproj);
    if (scale == NULL || scale[0] > 0.0) {
      if (scale == NULL)  for (i=0; i<nproj; i++) z[i] = x[proj[i] - 1];
      else {
	double invscale = 1.0 / scale[0];
	for (i=0; i<nproj; i++) z[i] = invscale * x[proj[i] - 1];
      }
    } else {
      // projection and aniso may not be given at the same time
      for (i=0; i<nproj; i++)
	z[i] = (x[proj[i] - 1] == 0 && scale[0] == 0) ? 0.0 : RF_INF;
    } 
    if (sign==NULL) COV(z, next, v) else LOGCOV(z, next, v, sign);
  } else if (Aniso != NULL) {
    int dim = Aniso->vdim[0];
    ALLOC_DOLLAR(z, dim);
    FCTN(x, Aniso, z);
    if (sign == NULL) COV(z, next, v) else LOGCOV(z, next, v, sign);
  } else if (aniso==NULL && (scale == NULL || scale[0] == 1.0)) {
    if (sign==NULL) COV(x, next, v) else LOGCOV(x, next, v, sign);
  } else {
    int xdimown = cov->xdimown;
    double *xz;
    ALLOC_DOLLAR(z, xdimown);
    if (aniso!=NULL) {
      xA(x, aniso, cov->nrow[DANISO], cov->ncol[DANISO], z);
      xz = z;
    } else xz = x;    
    if (scale != NULL) {
      if (scale[0] > 0.0) {
	double invscale = 1.0 / scale[0];
	for (i=0; i < xdimown; i++) z[i] = invscale * xz[i];
      } else {
	for (i=0; i < xdimown; i++)
	  z[i] = (xz[i] == 0.0 && scale[0] == 0.0) ? 0.0 : RF_INF;
      }
    }
    if (sign==NULL) COV(z, next, v) else LOGCOV(z, next, v, sign);
  }

  if (sign==NULL) {
    for (i=0; i<vdimSq; i++) v[i] *= var; 
  } else {
    double logvar = log(var);
    for (i=0; i<vdimSq; i++) v[i] += logvar; 
  }

}

void Snonstat(double *x, double *y, cov_model *cov, double *v){
  logSnonstat(x, y, cov, v, NULL);
}

void logSnonstat(double *x, double *y, cov_model *cov, double *v, 
		 double *sign){
  cov_model 
    *next = cov->sub[DOLLAR_SUB],
    *Aniso = cov->kappasub[DAUSER];
  double *z1, *z2,
    var = P0(DVAR),
    *scale =P(DSCALE),
    *aniso=P(DANISO);
  int i,
    nproj = cov->nrow[DPROJ],
    vdim = cov->vdim[0],
    vdimSq = vdim * vdim;

  //  PMI(cov, "$$$$");
  
  if (nproj > 0) {
    int *proj = PINT(DPROJ);
    ALLOC_DOLLARY(Z1, Z2, nproj);
    z1 = Z1;
    z2 = Z2;
    if (scale==NULL || scale[0] > 0.0) {
      double invscale = scale==NULL ? 1.0 :  1.0 / scale[0];
      for (i=0; i<nproj; i++) {
	z1[i] = invscale * x[proj[i] - 1];
	z2[i] = invscale * y[proj[i] - 1];	
      }
    } else {
      double s = scale[0]; // kann auch negativ sein ...
      for (i=0; i<nproj; i++) {
	z1[i] = (x[proj[i] - 1] == 0.0 && s == 0.0) ? 0.0 : RF_INF;
	z2[i] = (y[proj[i] - 1] == 0.0 && s == 0.0) ? 0.0 : RF_INF;
      }
    }
  } else if (Aniso != NULL) {
    int dim = Aniso->vdim[0];
    ALLOC_DOLLARY(Z1, Z2, dim);
    z1 = Z1;
    z2 = Z2;
    FCTN(x, Aniso, z1);
    FCTN(y, Aniso, z2);
    if (sign == NULL) NONSTATCOV(z1, z2, next, v)
      else LOGNONSTATCOV(z1, z2, next, v, sign);
  } else if (aniso==NULL && (scale==NULL || scale[0] == 1.0)) {
    z1 = x;
    z2 = y;
  } else {
    int xdimown = cov->xdimown;
    double *xz1, *xz2;
    ALLOC_DOLLARY(Z1, Z2, xdimown);
    z1 = Z1;
    z2 = Z2;
    if (aniso != NULL) {
      xA( x, y, aniso,cov->nrow[DANISO], cov->ncol[DANISO], z1, z2);
      xz1 = z1;
      xz2 = z2;
    } else {
      xz1 = x;
      xz2 = y;
    }
    if (scale != NULL) {
      double s = scale[0];
      if (s > 0.0) {
	double invscale = 1.0 / s;
	for (i=0; i<xdimown; i++) {
	  z1[i] = invscale * xz1[i];
	  z2[i] = invscale * xz2[i];
	}
      } else {
	for (i=0; i<nproj; i++) {
	  z1[i] = (xz1[i] == 0.0 && s == 0.0) ? 0.0 : RF_INF;
	  z2[i] = (xz2[i] == 0.0 && s == 0.0) ? 0.0 : RF_INF;
	}
      }
    }
  }
  
  if (sign == NULL) {
    NONSTATCOV(z1, z2, next, v);
    for (i=0; i<vdimSq; i++) v[i] *= var; 
  } else {
    double logvar = log(var);
    LOGNONSTATCOV(z1, z2, next, v, sign);
    for (i=0; i<vdimSq; i++) v[i] += logvar; 
  }
}


void covmatrixS(cov_model *cov, double *v, int *nonzeros) {
  location_type *loc = Loc(cov);	      
  cov_model *next = cov->sub[DOLLAR_SUB];
  location_type *locnext = Loc(next);
  assert(locnext != NULL);
  int i, err, tot, totSq,
    dim = loc->timespacedim,
    vdim = cov->vdim[0];
  double var = P0(DVAR);
  assert(dim == cov->tsdim);
    
  if ((err = alloc_cov(cov, dim, vdim, vdim)) != NOERROR) 
    error("memory allocation error in 'covmatrixS'");

  if ((!PisNULL(DSCALE) && P0(DSCALE) != 1.0) || 
      !PisNULL(DANISO) || !PisNULL(DPROJ) || 
      cov->kappasub[DSCALE] != NULL ||
      cov->kappasub[DAUSER] != NULL ||
       cov->kappasub[DPROJ] != NULL
      ) {
    CovarianceMatrix(cov, v, nonzeros); 
    return;
  }

  if (next->xdimprev != next->xdimown) {
    //printf("%d %d\n",  !PisNULL(DANISO), DANISO);
    //    PMI(cov, -1);
    BUG; // fuehrt zum richtigen Resultat, sollte aber nicht
    // vorkommen!
    CovarianceMatrix(cov, v, nonzeros); 
    return;
  }

  int next_gatter = next->gatternr,
    next_xdimgatter = next->xdimgatter,
    next_xdim = next->xdimprev;
 
  next->gatternr = cov->gatternr;
  next->xdimprev = cov->xdimprev;
  next->xdimgatter = cov->xdimgatter;
  CovList[next->nr].covmatrix(next, v, nonzeros);//hier wird uU next->totalpoints gesetzt
  next->gatternr = next_gatter;
  next->xdimgatter = next_xdimgatter;
  next->xdimprev = next_xdim;
  if (var == 1.0) return;

  // PMI(cov, "covmatrix S");
  // if (locnext==NULL) loc_set(cov, locnext->totalpoints);
  tot = cov->vdim[0] * locnext->totalpoints;
  totSq = tot * tot;
  for (i=0; i<totSq; v[i++] *= var);
}

char iscovmatrixS(cov_model *cov) {
  cov_model *sub = cov->sub[DOLLAR_SUB];
  return (int) ((PisNULL(DSCALE) || P0(DSCALE) ==1.0) &&	
		PARAMisNULL(cov, DAUSER) &&
		PARAMisNULL(cov, DPROJ) &&
		PARAMisNULL(cov, DANISO)) * CovList[sub->nr].is_covariance(sub);
}

void DS(double *x, cov_model *cov, double *v){
  cov_model *next = cov->sub[DOLLAR_SUB];
  assert(cov->kappasub[DAUSER] == NULL);
  int i,
    vdim = cov->vdim[0],
    vdimSq = vdim * vdim,
    nproj = cov->nrow[DPROJ];
  double y[2], varSc,
    *scale =P(DSCALE),
    *aniso=P(DANISO),
    spinvscale = 1.0;

  if (aniso != NULL) {
    spinvscale *= aniso[0];
    // was passiert, wenn es aniso nicht vom TypeIso ist ??
  }
  if (scale != NULL) spinvscale /= scale[0];
  varSc = P0(DVAR) * spinvscale;

  if (nproj == 0) {
    y[0] = x[0] * spinvscale; 
    y[1] = (cov->isoown==ISOTROPIC || cov->ncol[DANISO]==1) ? 0.0
      : x[1] * aniso[3]; // temporal; temporal scale
  } else {
    BUG;
    // for (i=0; i<nproj; i++) {
    //   y[i] = spinvscale * x[proj[i] - 1];
  }

  Abl1(y, next, v);
  for (i=0; i<vdimSq; i++) v[i] *= varSc; 
}

void DDS(double *x, cov_model *cov, double *v){
  cov_model *next = cov->sub[DOLLAR_SUB];
  assert(cov->kappasub[DAUSER] == NULL);
  int i,
    vdim = cov->vdim[0],
    vdimSq = vdim * vdim,
    nproj = cov->nrow[DPROJ],
    *proj = PINT(DPROJ);
  double y[2], varScSq,
    *scale =P(DSCALE),
    *aniso=P(DANISO),
    spinvscale = 1.0;
  
  if (aniso != NULL) {
    spinvscale *= aniso[0];
    // was passiert, wenn es aniso nicht vom TypeIso ist ??
  }
  if (scale != NULL) spinvscale /= scale[0];
  varScSq = P0(DVAR) * spinvscale * spinvscale;
  
  if (nproj == 0) {
    y[0] = x[0] * spinvscale;
    y[1] = (cov->isoown==ISOTROPIC || cov->ncol[DANISO]==1) ? 0.0
      : x[1] * aniso[3]; // temporal; temporal scale
  } else {
    BUG;
    for (i=0; i<nproj; i++) {
      y[0] += x[proj[i] - 1] * x[proj[i] - 1];
    }
    y[0] = sqrt(y[0]) * spinvscale;
  }
  Abl2(y, next, v);
  for (i=0; i<vdimSq; i++) v[i] *= varScSq; 
}


void D3S(double *x, cov_model *cov, double *v){
  cov_model *next = cov->sub[DOLLAR_SUB];
  assert(cov->kappasub[DAUSER] == NULL);
  int i,
    vdim = cov->vdim[0],
    vdimSq = vdim * vdim,
    nproj = cov->nrow[DPROJ],
    *proj = PINT(DPROJ);
  double y[2], varScS3,
    *scale =P(DSCALE),
    *aniso=P(DANISO),
    spinvscale = 1.0;

  if (aniso != NULL) {
    spinvscale *= aniso[0];
    // was passiert, wenn es aniso nicht vom TypeIso ist ??
  }
  if (scale != NULL) spinvscale /= scale[0];
  varScS3 = P0(DVAR) * spinvscale * spinvscale * spinvscale;
  
  if (nproj == 0) {
    y[0] = x[0] * spinvscale;
    y[1] = (cov->isoown==ISOTROPIC || cov->ncol[DANISO]==1) ? 0.0
      : x[1] * aniso[3]; // temporal; temporal scale
  } else {
    BUG;
    for (i=0; i<nproj; i++) {
      y[0] += x[proj[i] - 1] * x[proj[i] - 1];
    }
    y[0] = sqrt(y[0]) * spinvscale;
  }
  Abl3(y, next, v);
  for (i=0; i<vdimSq; i++) v[i] *= varScS3; 
}

void D4S(double *x, cov_model *cov, double *v){
  cov_model *next = cov->sub[DOLLAR_SUB];
  assert(cov->kappasub[DAUSER] == NULL);
  int i,
    vdim = cov->vdim[0],
    vdimSq = vdim * vdim,
    nproj = cov->nrow[DPROJ],
    *proj = PINT(DPROJ);
  double y[2], varScS4,
    *scale =P(DSCALE),
    *aniso=P(DANISO),
    spinvscale = 1.0;

  if (aniso != NULL) {
    spinvscale *= aniso[0];
    // was passiert, wenn es aniso nicht vom TypeIso ist ??
  }
  if (scale != NULL) spinvscale /= scale[0];
  varScS4 = spinvscale * spinvscale;
  varScS4 *= varScS4 * P0(DVAR);
  if (nproj == 0) {
    y[0] = x[0] * spinvscale;
    y[1] = (cov->isoown==ISOTROPIC || cov->ncol[DANISO]==1) ? 0.0
      : x[1] * aniso[3]; // temporal; temporal scale
  } else {
    BUG;
    for (i=0; i<nproj; i++) {
      y[0] += x[proj[i] - 1] * x[proj[i] - 1];
    }
    y[0] = sqrt(y[0]) * spinvscale;
  }
  Abl3(y, next, v);
  for (i=0; i<vdimSq; i++) v[i] *= varScS4; 
}


void nablahessS(double *x, cov_model *cov, double *v, bool nabla){
  cov_model *next = cov->sub[DOLLAR_SUB],
    *Aniso = cov->kappasub[DAUSER];
  if (Aniso != NULL) BUG;
  int i, endfor,
    dim = cov->nrow[DANISO],// == ncol == x d i m ?
    xdimown = cov->xdimown,
    nproj = cov->nrow[DPROJ];
  double *xy, *vw,
    *scale =P(DSCALE),
    *aniso=P(DANISO),
    var = P0(DVAR);
  if (nproj != 0) BUG;
  if (dim != xdimown) BUG;
	  
  
  if (aniso != NULL) {  
    ALLOC_DOLLAR(z, xdimown);
    ALLOC_DOLLAR2(w, xdimown);
    xA(x, aniso, xdimown, xdimown, z);
    xy = z;
    vw = w;
  } else {
    xy = x;
    vw = v;
  }

  if (scale != NULL) {
    ALLOC_DOLLAR3(y, xdimown);
    assert(scale[0] > 0.0);
    double spinvscale = 1.0 / scale[0];
    var *= spinvscale;
    if (!nabla) var *= spinvscale; // gives spinvscale^2 in case of hess
    for (i=0; i<xdimown; i++) y[i] = xy[i] * spinvscale;
    xy = y;
  }

  endfor = xdimown;
  if (nabla) {
    NABLA(xy, next, vw);
  } else {
    HESSE(xy, next, vw);
    endfor *= xdimown;
  }
     
  if (aniso != NULL) {  
    if (nabla) Ax(aniso, vw, xdimown, xdimown, v); 
    else {
      XCXt(aniso, vw, v, xdimown, xdimown); //todo:?reprogramm XCXt with alloc here ?
    }
  }

  for (i=0; i<endfor; i++) v[i] *= var; 
}

void nablaS(double *x, cov_model *cov, double *v){
  nablahessS(x, cov, v, true);
}
void hessS(double *x, cov_model *cov, double *v){
  nablahessS(x, cov, v, false);
}


 

void inverseS(double *x, cov_model *cov, double *v) {
  cov_model *next = cov->sub[DOLLAR_SUB];
  int i,
    idx[2] = {DAUSER, DPROJ};
  double 
    scale;
  
  for (i=0; i<2; i++) {
    if (cov->kappasub[idx[i]] != NULL) {
      char Msg[100];
      sprintf(Msg, "inverse can only be calculated if '%s' is not an arbitrary function",
  	      KNAME(idx[i])); 
      ERR(Msg);
    }
  }
  if (cov->kappasub[DSCALE] != NULL) {
    double left;
    NONSTATINVERSE(ZERO, cov->kappasub[DSCALE], &left, &scale);
    if (left < 0.0) ERR("scale not defined to be non-negative.");
  } else scale = PisNULL(DSCALE) ? 1.0 : P0(DSCALE);
 
  int
    dim = cov->xdimown,
    nproj = cov->nrow[DPROJ];
  //    *proj = (int *)P(DPROJ];
  double y, 
    s = 1.0,
    *aniso=P(DANISO),
    var = P0(DVAR);

  if (dim != 1) BUG;

  if (aniso != NULL) {
    if (isMiso(Type(aniso, cov->nrow[DANISO], cov->ncol[DANISO]))) 
      s /= aniso[0];
    else NotProgrammedYet(""); // to do
  }
  s *= scale;  
  if (nproj == 0) {
    y= *x / var; // inversion, so variance becomes scale
  } else {
    BUG;  //ERR("nproj is not allowed in invS");
  }
  
  if (CovList[next->nr].inverse == ErrCov) BUG;
  //PMI(next);
  INVERSE(&y, next, v);
 
  for (i=0; i<dim; i++) v[i] *= s; //!
}


void nonstatinverseS(double *x, cov_model *cov, double *left, double*right,
		     bool log){
  cov_model
    *next = cov->sub[DOLLAR_SUB],
    *Aniso = cov->kappasub[DAUSER];
   int i,
    dim = cov->tsdim,
    nproj = cov->nrow[DPROJ];
  //    *proj = (int *)P(DPROJ];
  double y, 
    s = 1.0,
    *scale =P(DSCALE),
    *aniso=P(DANISO),
    var = P0(DVAR);
 
 

  if (nproj == 0) {
    y= *x / var; // inversion, so variance becomes scale
  } else {
    BUG;  //ERR("nproj is not allowed in invS");
  }
    
  if (CovList[next->nr].nonstat_inverse == ErrInverseNonstat) BUG;
  if (log) {
    NONSTATLOGINVERSE(&y, next, left, right);
  } else {
    NONSTATINVERSE(&y, next, left, right);
  }

 
  if (aniso != NULL) {
    if (isMiso(Type(aniso, cov->nrow[DANISO], cov->ncol[DANISO]))) s/=aniso[0];
    else {
      dollar_storage *S = cov->Sdollar;
      int  
	ncol = cov->ncol[DANISO],
	nrow = cov->nrow[DANISO],
	ncnr = ncol * nrow,
	bytes = ncnr * sizeof(double),
	size = ncol * sizeof(double);
      bool redo;
      if (ncol != nrow) BUG;
      if ((redo = S->save_aniso == NULL)) {
	S->save_aniso = (double *) MALLOC(bytes);
	S->inv_aniso = (double *) MALLOC(bytes);
      }
      ALLOC_DOLLAR4(LR, ncol);
      double *save = S->save_aniso,
	*inv = S->inv_aniso;
      if (!redo) {
	for (i=0; i<ncnr; i++) if ((redo = save[i] != P(DANISO)[i])) break;
      }
      if (redo) {
	MEMCOPY(save, P(DANISO), bytes);
	MEMCOPY(inv, P(DANISO), bytes);
	if (invertMatrix(inv, nrow) != NOERROR)
	  error("inversion of anisotropy matrix failed");
      }
      
      MEMCOPY(LR, right, size);
      xA(LR, inv, nrow, ncol, right);
      
      MEMCOPY(LR, left, size);
      xA(LR, inv, nrow, ncol, left);      
    }    
  }

  if (Aniso != NULL) {
    if (aniso != NULL) BUG;

    //  PMI(Aniso);

    if (CovList[Aniso->nr].inverse == ErrInverse) 
      error("inverse of anisotropy matrix function unknown");
    int 
      nrow = Aniso->vdim[0],
      ncol = Aniso->vdim[1],
      size = nrow * sizeof(double);
    //      printf("ncol %d %d\n", ncol, nrow);
    if (cov->xdimown != ncol || ncol != 1)
      error("anisotropy function not of appropriate form");
    ALLOC_DOLLAR4(LR, nrow);
    
    MEMCOPY(LR, right, size);
    INVERSE(LR, Aniso, right);
    

    //    printf("right %f->%f  %f->%f", LR[0], right[0], LR[1], right[1]);
      
    MEMCOPY(LR, left, size);
    INVERSE(LR, Aniso, left);
  }

  if (scale != NULL) s *= scale[0];  
  if (s != 1.0) {
    for (i=0; i<dim; i++) {
      left[i] *= s; //!
      right[i] *= s;
    //      printf("i=%d l=%f %f \n", i, left[i], right[i]);      
    }
  }

  //APMI(cov->calling);

}

void nonstatinverseS(double *x, cov_model *cov, double *left, double*right) {
  nonstatinverseS(x, cov, left, right, false);
}

void nonstat_loginverseS(double *x, cov_model *cov, double *left, double*right){
 nonstatinverseS(x, cov, left, right, true);
}

void coinitS(cov_model *cov, localinfotype *li) {
  cov_model *next = cov->sub[DOLLAR_SUB];
  if ( CovList[next->nr].coinit == NULL)
    ERR("# cannot find coinit -- please inform author");
  CovList[next->nr].coinit(next, li);
}
void ieinitS(cov_model *cov, localinfotype *li) {
  cov_model *next = cov->sub[DOLLAR_SUB];
  
  if ( CovList[next->nr].ieinit == NULL)
    ERR("# cannot find ieinit -- please inform author");
  CovList[next->nr].ieinit(next, li);
}

void tbm2S(double *x, cov_model *cov, double *v){
  cov_model *next = cov->sub[DOLLAR_SUB];
  cov_fct *C = CovList + next->nr; // kein gatternr, da isotrop
  double y[2],  *xy,
    *scale =P(DSCALE),
    *aniso = P(DANISO);
  assert(cov->kappasub[DAUSER] == NULL);
  
  assert(cov->nrow[DPROJ] == 0);
  if (aniso!=NULL) {
    if (cov->ncol[DANISO]==2) {  // turning layers
      y[0] = x[0] * aniso[0]; // spatial 
      y[1] = x[1] * aniso[3]; // temporal
      assert(aniso[1] == 0.0 && aniso[2] == 0.0);
      if (y[0] == 0.0) C->cov(y, next, v); 
      else C->tbm2(y, next, v);
    } else {
      assert(cov->ncol[DANISO]==1);
      if (cov->nrow[DANISO] == 1) { // turning bands
	y[0] = x[0] * aniso[0]; // purely spatial 
	C->tbm2(y, next, v);
      } else { // turning layers, dimension reduction
	if (P0(DANISO) == 0.0) {
	  y[0] = x[1] * aniso[1]; // temporal 
	  C->cov(y, next, v); 
	} else {
	  y[0] = x[0] * aniso[0]; // temporal 
	  C->tbm2(y, next, v);
	}
      }
    }
    xy = y;
  } else xy = x;

  if (scale != NULL) {
    double s = scale[0];
    if (s > 0) { 
      double invscale = 1.0 / s;
      if (cov->xdimown == 2){
	y[0] = xy[0] * invscale; // spatial 
	y[1] = xy[1] * invscale; // temporal
	if (y[0] == 0.0) C->cov(y, next, v); 
	else C->tbm2(y, next, v);
      } else {
	y[0] = xy[0] * invscale; // purely spatial 
	C->tbm2(y, next, v);
      }
    } else {
      y[0] = (s < 0.0 || xy[0] != 0.0) ? RF_INF : 0.0;
      if (cov->xdimown == 2)
	y[1] = (s < 0.0 || xy[1] != 0.0) ? RF_INF : 0.0;
      C->tbm2(y, next, v);
    }
  }
  *v *= P0(DVAR);
}


// TODO : Aniso=Matrix: direkte implementierung in S,
// sodass nicht ueber initS gegangen werden muss, sondern
//  e  < -  e^\top Aniso


int TaylorS(cov_model *cov) {
  cov_model 
    *next = cov->sub[DOLLAR_SUB],
    *sub = cov->key == NULL ? next : cov->key;
  int i;

  if (PisNULL(DPROJ) && PisNULL(DANISO)) {
    double scale = PisNULL(DSCALE) ? 1.0 : P0(DSCALE);
    cov->taylorN = sub->taylorN;  
    for (i=0; i<cov->taylorN; i++) {
      cov->taylor[i][TaylorPow] = sub->taylor[i][TaylorPow];
      cov->taylor[i][TaylorConst] = sub->taylor[i][TaylorConst] *
	P0(DVAR) * pow(scale, -sub->taylor[i][TaylorPow]);   
    }

    cov->tailN = sub->tailN;  
    for (i=0; i<cov->tailN; i++) {
      cov->tail[i][TaylorPow] = sub->tail[i][TaylorPow];
      cov->tail[i][TaylorExpPow] = sub->tail[i][TaylorExpPow];
      cov->tail[i][TaylorConst] = sub->tail[i][TaylorConst] *
	P0(DVAR) * pow(scale, -sub->tail[i][TaylorPow]);   
      cov->tail[i][TaylorExpConst] = sub->tail[i][TaylorExpConst] *
	pow(scale, -sub->tail[i][TaylorExpPow]);
    }
  } else {
    cov->taylorN = cov->tailN = 0;
  }
  return NOERROR;
}

int checkS(cov_model *cov) {

  // hier kommt unerwartet  ein scale == nan rein ?!!
  cov_model 
    *next = cov->sub[DOLLAR_SUB],
    *Aniso = cov->kappasub[DAUSER],
    *sub = cov->key == NULL ? next : cov->key;
  int i, err,
    xdimown = cov->xdimown,
    xdimNeu = xdimown,
    *proj = PINT(DPROJ),
    nproj = cov->nrow[DPROJ];
  // bool skipchecks = GLOBAL.general.skipchecks;
  matrix_type type = TypeMany;

  
  assert(isAnyDollar(cov));
  if (!isDollarProc(cov)) cov->nr = DOLLAR; // wegen nr++ unten !
  
  // cov->q[1] not NULL then A has been given
  
  // if (cov->method == SpectralTBM && cov->xdimown != next->xdimprev)
  //  return ERRORDIM;
  
  if ((err = checkkappas(cov, false)) != NOERROR) {
    return err;
  }
  kdefault(cov, DVAR, 1.0);

  bool angle = Aniso!=NULL && CovList[Aniso->nr].check == checkAngle;
  if (angle && PisNULL(DANISO) && PisNULL(DAUSER) && cov->xdimown == cov->tsdim) {
    int dim  = cov->tsdim;
    if (!isCartesian(cov->isoown))  return ERRORANISO;
    if ((err = CHECK(Aniso, dim, dim, ShapeType, XONLY,
		    CARTESIAN_COORD, SUBMODEL_DEP, cov->role)) == NOERROR) {
      if (verysimple(Aniso)) {     
	 PALLOC(DAUSER, dim, dim);
 	 AngleMatrix(Aniso, P(DAUSER));
	 COV_DELETE(cov->kappasub + DAUSER);
	 Aniso = cov->kappasub[DAUSER] = NULL;
	 angle = false;
	 //PMI(cov, -1); 
      }
    }
  }


  if (cov->kappasub[DANISO] != NULL)
    SERR1("'%s' may not be a function.", KNAME(DANISO));
  if (!PisNULL(DAUSER)) {
    if (!isCartesian(cov->isoown))  return ERRORANISO;
    
    if (GLOBAL.internal.warn_Aniso) {
      PRINTF("NOTE! Starting with RandomFields 3.0, the use of '%s' is different from\nthe former '%s' insofar that '%s' is multiplied from the right by 'x' (i.e. Ax),\nwhereas '%s' had been multiplied from the left by 'x' (i.e. xA).\n", KNAME(DAUSER), KNAME(DANISO), KNAME(DANISO), KNAME(DAUSER));
    }
    GLOBAL.internal.warn_Aniso = false;
    // here for the first time
    if (!PisNULL(DANISO)) return ERRORANISO_T; 
    int k,
      lnrow = cov->nrow[DAUSER],
      lncol = cov->ncol[DAUSER];
    long j, 
      total = lncol * lnrow;
	
    double
      *pA = P(DAUSER); 
    PALLOC(DANISO, lncol, lnrow); // !! ACHTUNG col, row gekreuzt
    for (i=k=0; i<lnrow; i++) {
      for (j=i; j<total; j+=lnrow) P(DANISO)[k++] = pA[j];
    }
    PFREE(DAUSER);
  }
 

  if (Aniso != NULL) {

    if (!isDollarProc(cov) && !angle) cov->domown = KERNEL;

    if (!PisNULL(DANISO) || !PisNULL(DPROJ) || !PisNULL(DSCALE))
      SERR2("if '%s' is an arbitrary function, only '%s' may be given as additional parameter.", KNAME(DAUSER), KNAME(DVAR));
    if (cov->isoown != SYMMETRIC && cov->isoown != CARTESIAN_COORD) {
      //      printf("iso = %d\n", cov->isoown);
      //PMI(cov, -1);
      //      assert(false);
      return ERRORANISO;
    } //else  printf("OK iso = %d\n", cov->isoown);
 
    cov->full_derivs = cov->rese_derivs = 0;
    cov->loggiven = true;
    if ((err = CHECK(Aniso, cov->tsdim, cov->xdimown, ShapeType, XONLY,
		     CARTESIAN_COORD, SUBMODEL_DEP, cov->role)) != NOERROR) {
      return err;
    }

    if (Aniso->vdim[1] != 1)
      SERR4("'%s' returns a matrix of dimension %d x %d, but dimension %d x 1 is need.",
	    KNAME(DANISO), Aniso->vdim[0], Aniso->vdim[1], cov->xdimown);

    // PMI(cov);

    if (cov->key==NULL) {
      if ((err = CHECK(sub, Aniso->vdim[0], Aniso->vdim[0], cov->typus,  
		       cov->domown, 
		       cov->isoown, SUBMODEL_DEP, cov->role)) != NOERROR) {
	return err;
      }
    }
    cov->pref[Nugget] = cov->pref[RandomCoin] = cov->pref[Average] = 
      cov->pref[Hyperplane] = cov->pref[SpectralTBM] = cov->pref[TBM] = 
      PREF_NONE;

    sub->pref[CircEmbed] = sub->pref[CircEmbedCutoff] = 
      sub->pref[CircEmbedIntrinsic] = sub->pref[Sequential] = 
      sub->pref[Specific] = PREF_NONE;

 
  } else if (!PisNULL(DANISO)) { // aniso given
    int 
      nrow = cov->nrow[DANISO],
      ncol = cov->ncol[DANISO];
 
    if (nrow==0 || ncol==0) SERR("dimension of the matrix is 0");
    if (!PisNULL(DPROJ)) return ERRORANISO_T;
    if (xdimown < nrow) {
      if (PL >= PL_ERRORS) {LPRINT("xdim=%d != nrow=%d\n", xdimown, nrow);}
      SERR("#rows of anisotropy matrix does not match dim. of coordinates");
    }
    if (xdimown != cov->tsdim && nrow != ncol)
      SERR("non-quadratic anisotropy matrices only allowed if dimension of coordinates equals spatio-temporal dimension");
 
    type = Type(P(DANISO), nrow, ncol);
    
    cov->full_derivs = cov->rese_derivs = 0;
    cov->loggiven = true;


    switch (cov->isoown) {
    case ISOTROPIC :   
      if (cov->tsdim != 1) return ERRORANISO;
      cov->full_derivs = cov->rese_derivs = 2;
      break;
    case SPACEISOTROPIC :  
      cov->full_derivs =  cov->rese_derivs = isMdiag(type) ? 2 : 0;
      if (nrow != 2 || !isMdiag(type)) {
	//	printf("nrow %d %d %d\n", nrow, type, isMdiag(type));
	//	int h; for (h = 0; h<cov->nrow[DANISO] * cov->ncol[DANISO]; h++)
	//		 printf("h=%d %f %d\n", h, P(DANISO)[h], P(DANISO)[h]==0);
	SERR("spaceisotropy needs a 2x2 diagonal matrix");
      }
      break;      
    case ZEROSPACEISO : 
      if (!isMdiag(type)) return ERRORANISO;
      break;      
    case VECTORISOTROPIC :
      if (!isMiso(type)) return ERRORANISO; 
      break;
    case SYMMETRIC: 
      break;
    case PREVMODELI : BUG;      
    case CARTESIAN_COORD : case GNOMONIC_PROJ :  case ORTHOGRAPHIC_PROJ:
      if (!isProcess(cov->typus)) return ERRORANISO;
      break;
    default : 
      if (isCartesian(cov->isoown)) { BUG; }
      else return  ERRORANISO;
    }

    //  printf(" A hiere\n");
    
    if (!isMdiag(type)) 
      cov->pref[Nugget] = cov->pref[RandomCoin] = cov->pref[Average] = cov->pref[Hyperplane] = PREF_NONE;
    if (cov->isoown != SPACEISOTROPIC && !isMiso(type))
      cov->pref[SpectralTBM] = cov->pref[TBM] = PREF_NONE;
   
    //    PMI(cov);
    //printf("ncol=%d %d %s %s\n", ncol, isProcess(cov->typus),  NICK(cov), 
    //	   ISONAMES[cov->isoown]);
    if (cov->key==NULL) {
      if ((err = CHECK(next, ncol, ncol, cov->typus, cov->domown, 
		 ncol==1 && !isProcess(cov->typus) ? ISOTROPIC : cov->isoown, 
		 SUBMODEL_DEP, cov->role))
	  != NOERROR) {
	// assert(false); 
 	return err;
      }
      if (next->domown == cov->domown && next->isoown == cov->isoown &&
	  xdimown > 1) next->delflag = DEL_COV - 7;
    } else {
      if ((err = CHECK(cov->key, ncol, ncol, cov->typus, cov->domown, 
		       ncol==1 && 
		       !isProcess(cov->typus) ? ISOTROPIC : cov->isoown,
		       SUBMODEL_DEP, cov->role)) != NOERROR) return err;
    }
  
  } else { // PisNULL(DANISO) 
    
    int tsdim = cov->tsdim;
    if (nproj > 0) {
      if (cov->ncol[DPROJ] != 1) SERR("proj must be a vector");
      if (cov->xdimprev != cov->xdimown) return ERRORANISO;
      for (i=0; i<nproj; i++) {
	int idx = proj[i] - 1;
	if (idx >= xdimown)
	  SERR3("%d-th value of '%s' (%d) out of range", 
		i, KNAME(DPROJ), proj[i]);
      }
      xdimNeu = nproj;
      tsdim = nproj;

      switch (cov->isoown) {
      case ISOTROPIC : if (cov->tsdim != 1) return ERRORANISO;
	break;
      case SPACEISOTROPIC :  
	if (nproj != 2 || xdimown != 2 || tsdim != 2)
	  SERR("spaceisotropy needs a 2x2 diagonal matrix");
	break;      
      case ZEROSPACEISO : return ERRORANISO; // ginge z.T.; aber kompliziert
	break;      
      case VECTORISOTROPIC : return ERRORANISO; 
	break;
      case SYMMETRIC: case CARTESIAN_COORD: case GNOMONIC_PROJ :
      case ORTHOGRAPHIC_PROJ :
	break;
      case PREVMODELI : BUG;      
	break;
      default : 
	if (isCartesian(cov->isoown)) {BUG;}
	else return  ERRORANISO;  // todo
      }
    }

    if (cov->key==NULL) {
      //PMI(next);
      if ((err = CHECK(next, tsdim, xdimNeu, cov->typus, cov->domown,
		       cov->isoown, 
		       cov->vdim[0], // SUBMODEL_DEP; geaendert 20.7.14
		       cov->role)) != NOERROR) {
	return err;
      }

      if (next->domown == cov->domown &&
	  next->isoown == cov->isoown) // darf nicht nach C-HECK als allgemeine Regel ruebergezogen werden, u.a. nicht wegen stat/nicht-stat wechsel !!
	// hier geht es, da domown und isoown nur durchgegeben werden und die Werte      // bereits ein Schritt weiter oben richtig/minimal gesetzt werden.
	next->delflag = DEL_COV - 8;
    } else {
      if ((err = CHECK(cov->key, tsdim, xdimNeu, cov->typus,
		       cov->domown, cov->isoown,
		       SUBMODEL_DEP, cov->role)) != NOERROR) return err;
    }

  } // end no aniso
 
  //  printf("tsdim = %d\n", cov->tsdim);

  if (( err = checkkappas(cov, false)) != NOERROR) {
    return err;
  }


  //     printf("here %d\n", cov->maxdim);
  
  setbackward(cov, sub);
  if ((Aniso != NULL || !PisNULL(DANISO) || !PisNULL(DPROJ)) && 
      cov->maxdim < cov->xdimown) cov->maxdim = cov->xdimown;
 	
  if (!isAnyIsotropic(cov->isoown) && !isDollarProc(cov)) { // multivariate kann auch xdimNeu == 1 problematisch sein
    cov->nr++;
  }
  
  if (xdimNeu > 1) {
    cov->pref[CircEmbedCutoff] = cov->pref[CircEmbedIntrinsic] = 0;
  }

  // printf("cov-Nr = %d\n", cov->nr);
  
  // 30.10.11 kommentiert:
  //  cov->pref[CircEmbedCutoff] = cov->pref[CircEmbedIntrinsic] = 
  //      cov->pref[TBM] = cov->pref[SpectralTBM] = 0;
  if ( (PisNULL(DANISO) || isMiso(type)) && PisNULL(DPROJ)) {
    //    print("logspeed %s %f\n", NICK(sub), sub->logspeed);
    cov->logspeed = sub->logspeed * P0(DVAR);
  }
  //////////////// 

  if (sub->pref[Nugget] > PREF_NONE) cov->pref[Specific] = 100;
  if (nproj == 0) cov->matrix_indep_of_x = sub->matrix_indep_of_x;

  if ((err = TaylorS(cov)) != NOERROR) return err;
  // printf("here2\n");
  
  if (cov->finiterange && isAnySpherical(cov->isoown)) {
    double range;
    INVERSE(ZERO, cov, &range);
    bool ok = (isSpherical(cov->isoown) && range <= M_PI) ||
      (isEarth(cov->isoown) && range <= 180);
    if (!ok)
      SERR2("model '%s' does not allow for required coordinate system '%s'",
	    NICK(next), COORD_SYS_NAMES[cov->isoown]);
  }
  
  DOLLAR_STORAGE;
  
  if (isProcess(cov->typus)) {
    MEMCOPY(cov->pref, PREF_NOTHING, sizeof(pref_shorttype)); 
  }
 
  if (GLOBAL.coords.coord_system == earth && 
      is_all(isCartesian, CovList + next->nr) &&
      GLOBAL.internal.warn_scale &&
      P0(DSCALE) < (strcmp(GLOBAL.coords.newunits[0], "km") == 0 ? 10 : 6.3)) {
    char msg[300];

    //  printf("%d %d\n", P0(DSCALE), strcmp(GLOBAL.general.newunits[0], "km") == 0 ? 100 : 63); assert(false);

    sprintf(msg, "value of scale parameter equals '%4.2f',\nwhich is less than 100, although models defined on R^3 are used in the\ncontext of earth coordinates where larger scales are expected.\n(This warning appears only ones per session.)\n", P0(DSCALE)); 
    GLOBAL.internal.warn_scale = false;
    warning(msg);
  }
  //  PMI(cov, 1);
  //  printf("here XX\n");
  //PMI(cov->calling);
   return NOERROR;
}



void rangeS(cov_model *cov, range_type* range){
  int i;

  range->min[DVAR] = 0.0;
  range->max[DVAR] = RF_INF;
  range->pmin[DVAR] = 0.0;
  range->pmax[DVAR] = 100000;
  range->openmin[DVAR] = false;
  range->openmax[DVAR] = true;

  range->min[DSCALE] = 0.0;
  range->max[DSCALE] = RF_INF;
  range->pmin[DSCALE] = 0.0001;
  range->pmax[DSCALE] = 10000;
  range->openmin[DSCALE] = true;
  range->openmax[DSCALE] = true;

  for (i=DANISO; i<= DAUSER; i++) {
    range->min[i] = RF_NEGINF;
    range->max[i] = RF_INF;
    range->pmin[i] = -10000;
    range->pmax[i] = 10000;
    range->openmin[i] = true;
    range->openmax[i] = true;
  }
  
  range->min[DPROJ] = 1;
  range->max[DPROJ] = cov->tsdim;
  range->pmin[DPROJ] = 1;
  range->pmax[DPROJ] =  cov->tsdim;
  range->openmin[DPROJ] = false;
  range->openmax[DPROJ] = false;
}


bool TypeS(Types required, cov_model *cov, int depth) {
  cov_model *sub = cov->key==NULL ? cov->sub[0] : cov->key;

  // printf("\n\n\nxxxxxxxxxxxx %d %s %s\n\n", TypeConsistency(required, sub, 0), TYPENAMES[required], NICK(sub));
  
  // if (required == ProcessType) crash(); //assert(false);

  return (required==TcfType || required==PosDefType || required==NegDefType
	  || required==ShapeType || required==TrendType 
	  || required==ProcessType || required==GaussMethodType)
    && TypeConsistency(required, sub, depth-1);
}


void spectralS(cov_model *cov, gen_storage *s, double *e) {
  cov_model *next = cov->sub[DOLLAR_SUB];
  int d,
    ncol = PisNULL(DANISO) ? cov->tsdim : cov->ncol[DANISO];
  double sube[MAXTBMSPDIM],
    *scale =P(DSCALE),
    invscale = 1.0;
  
  SPECTRAL(next, s, sube); // nicht gatternr

  // Reihenfolge nachfolgend extrem wichtig, da invscale auch bei aniso
  // verwendet wird


  if (scale != NULL) invscale /= scale[0];
  // print("sube %f %f %f %d %d\n", sube[0], sube[1], invscale, cov->tsdim, ncol);
  
  if (!PisNULL(DANISO)) {
    int 
      nrow = cov->nrow[DANISO];
    long j, k, m,
      total = ncol * nrow;
    double
      *A = P(DANISO); 
    for (d=0, k=0; d<nrow; d++, k++) {
      e[d] = 0.0;
      for (m=0, j=k; j<total; j+=nrow) {
	e[d] += sube[m++] * A[j] * invscale;
      }
    }
  } else { 
    for (d=0; d<ncol; d++) e[d] = sube[d] * invscale;
  }

}


void ScaleDollarToLoc(cov_model *to, cov_model *from,
		      int VARIABLE_IS_NOT_USED depth) {

  assert(!PARAMisNULL(to, LOC_SCALE));
  assert(isDollar(from));
  assert(!PARAMisNULL(from, DSCALE));
  PARAM(to, LOC_SCALE)[0] = PARAM0(from, DSCALE);
}

bool ScaleOnly(cov_model *cov) {
  return isDollar(cov) && 
    PisNULL(DPROJ) && cov->kappasub[DPROJ] == NULL &&
    PisNULL(DANISO) &&  cov->kappasub[DAUSER] == NULL &&
    (PisNULL(DVAR) || P0(DVAR)==1.0) && cov->kappasub[DVAR] == NULL;
}



int addScales(cov_model **newmodel, double anisoScale, cov_model *scale,
	      double Scale) {
  //  cov_model *calling = (*newmodel)->calling;

  if (anisoScale != 1.0) {
    addModel(newmodel, LOC);
    kdefault(*newmodel, LOC_SCALE, anisoScale);
  }

  if (scale != NULL) {
    if (!isRandom(scale)) SERR("unstationary scale not allowed yet");
    addModel(newmodel, LOC);
    addSetDistr(newmodel, scale->calling, ScaleDollarToLoc, true, MAXINT);

    //   assert((*newmodel)->Sset != NULL);
    // APMI(*newmodel);
    
  } else {
    if (Scale != 1.0) {
      addModel(newmodel, LOC);
      kdefault(*newmodel, LOC_SCALE, Scale);
    }
  }  
  return NOERROR;
}

int structS(cov_model *cov, cov_model **newmodel) {
  if (cov->role == ROLE_GAUSS && isProcess(cov->typus)) {
    cov->nr = DOLLAR_PROC;
    return structSproc(cov, newmodel); // kein S-TRUCT(...) !!
  } 

  cov_model *local = NULL,
    *dummy = NULL,
    *next = cov->sub[DOLLAR_SUB],
    *Aniso = cov->kappasub[DAUSER],
    *scale = cov->kappasub[DSCALE];
  int err = NOERROR; 
  bool generalAniso = Aniso != NULL && CovList[Aniso->nr].check != checkAngle;

   
  ASSERT_NEWMODEL_NOT_NULL;
 
  if (cov->kappasub[DVAR] != NULL) {
    GERR1("models including arbitrary functions for '%s' cannot be simulated yet", KNAME(DVAR));
  } 

  if (generalAniso)
    GERR1("complicated models including arbitrary functions for '%s' cannot be simulated yet", KNAME(DAUSER));

  //APMI(cov->calling->calling->calling);
  
  switch (cov->role) {
  case ROLE_POISSON_GAUSS : 
  case ROLE_SMITH : 
    if (!next->deterministic) GERR("random shapes not programmed yet");
    if (!PisNULL(DANISO)) GERR("anisotropy parameter not allowed yet");
    if (!PisNULL(DPROJ)) GERR("projections not allowed yet");
    if ((err = STRUCT(next, newmodel)) > NOERROR) return err;

    addModel(newmodel, DOLLAR);
    if (!PisNULL(DVAR)) kdefault(*newmodel, DVAR, P0(DVAR));
    if (!PisNULL(DSCALE)) kdefault(*newmodel, DSCALE, P0(DSCALE));
    break;

  case ROLE_MAXSTABLE : // eigentlich nur von RPSmith moeglich !
    if (!next->deterministic) GERR("random shapes not programmed yet");
    if (!PisNULL(DPROJ)) GERR("projections not allowed yet");
    // P(DVAR) hat keine Auswirkungen
    if (!PisNULL(DANISO)) {
      if (!isMonotone(next) || !isIsotropic(next->isoown) ||
	  PisNULL(DANISO) || cov->ncol[DANISO] != cov->nrow[DANISO])
	GERR("anisotropy parameter only allowed for simple models up to now");
    }

    //APMI(cov);
    
    assert(cov->calling != NULL); 
    double anisoScale;
  
    if (!PisNULL(DANISO)) {
      anisoScale = 1.0 / getMinimalAbsEigenValue(P(DANISO), cov->nrow[DANISO]);
      if (!R_FINITE(anisoScale)) 
	GERR("singular anisotropy matrices not allowed");
    } else anisoScale = 1.0;
  
    if (cov->calling->nr == SMITHPROC) {
      if ((err = STRUCT(next, newmodel)) == NOERROR && *newmodel != NULL) {
	assert(	(*newmodel)->calling == cov);
	//   PMI(shape); 
	//APMI(*newmodel);
	Types type = 
	  TypeConsistency(PointShapeType, *newmodel, 0) ? PointShapeType :
	  TypeConsistency(RandomType, *newmodel, 0) ? RandomType :
	  TypeConsistency(ShapeType, *newmodel, 0) ? ShapeType :
	  OtherType;
	double Scale =  PisNULL(DSCALE) ? 1.0 : P0(DSCALE);
	
	if (type == RandomType) { // random locations given;
	  // so, it must be of pgs type (or standard):
	  
	  if ((err = CHECK_R(*newmodel, cov->tsdim)) != NOERROR) {
	    goto ErrorHandling;
	  }
	  dummy = *newmodel;

	  if ((err=addScales(&dummy, anisoScale, scale, Scale))!=NOERROR){
	    goto ErrorHandling;
	  }	  

	  *newmodel = NULL;     
	  assert(cov->vdim[0] == cov->vdim[1]);

	  
	  //APMI(dummy);

	  if ((err = addPointShape(newmodel, cov, dummy, cov->tsdim, 
				   cov->vdim[0])) != NOERROR) {
	    goto ErrorHandling; 
	  }

	  ASSERT_NEWMODEL_NOT_NULL;
	  (*newmodel)->calling = cov;
	  // APMI(cov);
	} else { // type not RandomType
	  if (type == PointShapeType && 
	      (err = addScales((*newmodel)->sub + PGS_LOC, anisoScale, scale, 
			       Scale)) != NOERROR) goto ErrorHandling;
	  if ((err = CHECK(*newmodel, cov->tsdim, cov->xdimprev, type, 
			   cov->domprev, cov->isoprev, cov->vdim, 
			   ROLE_MAXSTABLE))
	      != NOERROR) {
	    goto ErrorHandling;
	  }
	  if (type==PointShapeType) {	
	    if ((err = PointShapeLocations(*newmodel, cov)) != NOERROR) 
	      goto ErrorHandling;
	  } else if (type == ShapeType) {
	    dummy = *newmodel;
	    *newmodel = NULL;
	    // suche nach geeigneten locationen
	    if ((err = addScales(&local, anisoScale, scale, Scale))!=NOERROR)
	      goto ErrorHandling;
	    if ((err = addPointShape(newmodel, dummy, NULL, local,
				     cov->tsdim, cov->vdim[0]))
		!= NOERROR) 
	      goto ErrorHandling; 
	    //printf("e ddd\n");
	  } else BUG;
	} // ! randomtype
      } else { // S TRUCT does not return anything
	int err2;
	//assert(false);
	// XERR(err);
	//		APMI(*newmodel);
	if ((err2 = addPointShape(newmodel, cov, NULL, cov->tsdim,
				  cov->vdim[0]))
	    != NOERROR) {
	  if (err == NOERROR) err = err2;
	  goto ErrorHandling; 
	}
	err = NOERROR;
      }

    } else { // not from RPsmith
      BUG;
      // 
     if ((err = STRUCT(next, newmodel)) > NOERROR) return err;
    }
    
    //  PMI(cov);

    //  APMI(*newmodel);
   
    break;
  case ROLE_GAUSS :
    if (cov->key != NULL) COV_DELETE(&(cov->key));

    if (cov->prevloc->distances) 
      GERR("distances do not allow for more sophisticated simulation methods");
    
    if ((err = STRUCT(next, newmodel)) > NOERROR) return err;

    addModel(newmodel, DOLLAR);
    if (!PisNULL(DVAR)) kdefault(*newmodel, DVAR, P0(DVAR));
    if (!PisNULL(DSCALE) ) kdefault(*newmodel, DSCALE, P0(DSCALE));
    if (!next->deterministic) GERR("random shapes not programmed yet");
    
    break;
 
  default :
    //  PMI(cov, "structS");
    GERR2("%s : changes in scale/variance not programmed yet for '%s'", 
	  NICK(cov), ROLENAMES[cov->role]);      
    }
  
  ErrorHandling:
  
    if (dummy != NULL) COV_DELETE(&dummy);
    if (local != NULL) COV_DELETE(&local);

   
  return err;
}




int initS(cov_model *cov, gen_storage *s){
  // am liebsten wuerde ich hier die Koordinaten transformieren;
  // zu grosser Nachteil ist dass GetDiameter nach trafo 
  // grid def nicht mehr ausnutzen kann -- umgehbar?!
  cov_model *next = cov->sub[DOLLAR_SUB],
    *varM = cov->kappasub[DVAR],
    *scaleM = cov->kappasub[DSCALE],
    *anisoM = cov->kappasub[DAUSER],
    *projM = cov->kappasub[DPROJ];
  int 
    vdim = cov->vdim[0],
    nm = cov->mpp.moments,
    nmvdim = (nm + 1) * vdim,
    err = NOERROR;
  bool 
    angle = anisoM != NULL && CovList[anisoM->nr].check == checkAngle,
    maxstable = hasExactMaxStableRole(cov);// Realisationsweise 

  if (hasAnyShapeRole(cov)) { // Average  !! ohne maxstable selbst !!
    double
      var[MAXMPPVDIM],
      scale = PisNULL(DSCALE)  ? 1.0 : P0(DSCALE);
    int i,
      dim = cov->tsdim;
    

    if (!PisNULL(DPROJ) || !PisNULL(DANISO) || 
	projM!= NULL || (anisoM != NULL && (!angle || !anisoM->deterministic))){
       //      APMI(cov);

       SERR("(complicated) anisotropy ond projection not allowed yet in Poisson related models");
     }


    // Achtung I-NIT_RANDOM ueberschreibt mpp.* !!
    if (varM != NULL) {
      int nm_neu = nm == 0 && !maxstable ? 1 : nm;
      if ((err = INIT_RANDOM(varM, nm_neu, s, P(DVAR))) != NOERROR) return err; 
      
      int nmP1 = varM->mpp.moments + 1;
      for (i=0; i<vdim; i++) {
	int idx = i * nmP1;
	var[i] = maxstable ? P0(DVAR) : varM->mpp.mM[idx + 1];      
      }
    } else for (i=0; i<vdim; var[i++] = P0(DVAR));

    if (scaleM != NULL) {
      if (dim + nm < 1) SERR("found dimension <= 0");
      int dim_neu = maxstable ? nm : (dim + nm) < 1 ? 1 : dim + nm; 
    if ((err = INIT_RANDOM(scaleM, dim_neu, s, P(DSCALE)))
	!= NOERROR) return err;
      scale = maxstable ? P0(DSCALE) : scaleM->mpp.mM[1];      
    }
    if ((err = INIT(next, nm, s)) != NOERROR) return err;


    for (i=0; i < nmvdim; i++) {
      cov->mpp.mM[i] = next->mpp.mM[i]; 
      cov->mpp.mMplus[i] = next->mpp.mMplus[i]; 
    }

    if (varM != NULL && !maxstable) {
      for (i=0; i < nmvdim; i++) {
	cov->mpp.mM[i] *= varM->mpp.mM[i];
	cov->mpp.mMplus[i] *= varM->mpp.mMplus[i];
      }
    } else {
      int j, k;
      double pow_var;
      for (k=j=0; j<vdim; j++) { 
	pow_var = 1.0;
	for (i=0; i <= nm; i++, pow_var *= var[j], k++) {	
	  cov->mpp.mM[k] *= pow_var;
	  cov->mpp.mMplus[k] *= pow_var;
	}
      }
    }
    if (scaleM != NULL && !maxstable) {
      if (scaleM->mpp.moments < dim) SERR("moments can not be calculated.");
      int j, k,
	nmP1 = scaleM->mpp.moments + 1;
      for (k=j=0; j<vdim; j++) { 
	double pow_scale = scaleM->mpp.mM[dim + j * nmP1];
	for (i=0; i <= nm; i++, k++) {
	  cov->mpp.mM[k] *= pow_scale;
	  cov->mpp.mMplus[k] *= pow_scale;
	}
      }
    } else {
      double pow_scale = pow(scale, dim);
      for (i=0; i < nmvdim; i++) {
	cov->mpp.mM[i] *= pow_scale;
	cov->mpp.mMplus[i] *= pow_scale;
      }
    }

    if (!PisNULL(DANISO)) {      
      if (cov->nrow[DANISO] != cov->ncol[DANISO])
	SERR("only square anisotropic matrices allowed");
      double invdet = fabs(1.0 / getDet(P(DANISO), cov->nrow[DANISO]));
      if (!R_FINITE(invdet))
	SERR("determinant of the anisotropy matrix could not be calculated.");
      
      for (i=0; i < nmvdim; i++) {
	cov->mpp.mM[i] *= invdet;
	cov->mpp.mMplus[i] *= invdet;
      }

    }

    if (anisoM != NULL) {  
      double invdet;
      if (angle) {
	int 
	  ncol = anisoM->vdim[0],
	  nrow = anisoM->vdim[1];
	if (nrow != ncol) SERR("only square anisotropic matrices allowed");
	double 
	  *diag = PARAM(anisoM, ANGLE_DIAG);
	if (diag != NULL) {
	  invdet = 1.0;
	  for (i=0; i<ncol; i++) invdet /= diag[i];
	} else {
	  invdet = PARAM0(anisoM, ANGLE_RATIO);
	}
      } else {
	SERR("general anisotropy matrices not allowed yet");
      }
      
      invdet = fabs(invdet);      
      if (!R_FINITE(invdet))
	SERR("determinant of the anisotropy matrix could not be calculated.");
      
      for (i=0; i < nmvdim; i++) {
	cov->mpp.mM[i] *= invdet;
	cov->mpp.mMplus[i] *= invdet;
      }
    }

    
 
    if (R_FINITE(cov->mpp.unnormedmass)) {
      if (vdim > 1) BUG;
      cov->mpp.unnormedmass = next->mpp.unnormedmass * var[0];
    } else {
      for (i=0; i<vdim; i++)
	cov->mpp.maxheights[i] = next->mpp.maxheights[i] * var[i];
    }
 
    //    printf("inis %f %f %d %f\n", cov->mpp.maxheights[0] , next->mpp.maxheights[0] , varM == NULL, varM == NULL ? P0(DVAR) : 1.0);
 
    //    cov->mpp.refradius *= scale;
    //cov->mpp.refsd *= scale;
  }


  else if (cov->role == ROLE_GAUSS) {
    cov_model 
      *key = cov->key,
      *sub = key == NULL ? next : key;
    assert(sub != NULL);

    assert(key == NULL || ({PMI(cov);false;}));//
      
   
    if ((err=INIT(sub, 0, s)) != NOERROR) return err;
   
  }

  else if (cov->role == ROLE_BASE) {
    cov_model 
      *key = cov->key,
      *sub = key == NULL ? next : key;
    assert(sub != NULL);

    assert(key == NULL || ({PMI(cov);false;}));//
      
   
    if ((err=INIT(sub, 0, s)) != NOERROR) return err;
    
  } else SERR("initiation of scale and/or variance failed");

 if ((err = TaylorS(cov)) != NOERROR) return err;

  //  APMI(cov);

 return NOERROR;

}


void doS(cov_model *cov, gen_storage *s){
   cov_model
     *varM = cov->kappasub[DVAR],
     *scaleM = cov->kappasub[DSCALE];
   int i,
     vdim = cov->vdim[0];

    if (varM != NULL && !varM->deterministic && !isRandom(varM)) {
      assert(!PisNULL(DVAR));
      DO(varM, s);
    }
    
    if (scaleM != NULL && !scaleM->deterministic && !isRandom(scaleM)) {
      assert(!PisNULL(DSCALE));
      DO(scaleM, s);
    }

    //    PMI(cov);


  if (hasAnyShapeRole(cov)) {
    cov_model *next = cov->sub[DOLLAR_SUB];
         
    DO(next, s);// nicht gatternr
    for (i=0; i<vdim; i++)
      cov->mpp.maxheights[i] = next->mpp.maxheights[i] * P0(DVAR);
    return;
  }  
  
  else if (cov->role == ROLE_GAUSS) {    
    double 
      *res = cov->rf,
      sd = sqrt(P0(DVAR));
    int 
      totalpoints = Gettotalpoints(cov);
    assert(res != NULL);
    if (cov->key == NULL) BUG;

    DO(cov->key, s); 

    //    PMI(cov);
    //    printf("totalpoints");

    if (sd != 1.0) for (i=0; i<totalpoints; i++) res[i] *= (res_type) sd;
    return;
  } 

  //PMI(cov->calling->calling);
  //crash();

  BUG;
}



int checkplusmal(cov_model *cov) {
  cov_model *sub;
  int i, j, err, dim, xdim, role;
  
  assert(cov->Splus == NULL);  
  dim = cov->tsdim;
  xdim = cov->xdimown;
  role = cov->role;

  //  PMI(cov, -1);
  //  assert(cov->typus == PosDefType);
  
  for (i=0; i<cov->nsub; i++) { 

    sub = cov->sub[i];
     
    if (sub == NULL) 
      SERR("+ or *: named submodels are not given in a sequence!");

    Types type = cov->typus;
    if (CovList[cov->nr].check == checkmal && type==NegDefType) type=PosDefType;
    domain_type dom = cov->domown;
    isotropy_type iso = cov->isoown;

    //    printf("start\n");
    err = ERRORTYPECONSISTENCY;
    for (j=0; j<=1; j++) { // nur trend als abweichender typus erlaubt
      if (TypeConsistency(type, sub, 0) &&
	  (err = CHECK(sub, dim, xdim, type, dom, iso, 
		       i == 0 ? SUBMODEL_DEP : cov->vdim[0], role))
	  == NOERROR) break;
      type = TrendType;
      dom = XONLY;
      iso = CARTESIAN_COORD;
    }
    
    if (err != NOERROR) {
      return err;
    }

    if (cov->typus == sub->typus) {
      setbackward(cov, sub);
    } else {  
      updatepref(cov, sub);
      cov->tbm2num |= sub->tbm2num;
      if (CovList[cov->nr].vdim == SUBMODEL_DEP && 
	  (sub==cov->sub[0] || sub==cov->key)) { // strange todo
	cov->vdim[0] = sub->vdim[0];
	cov->vdim[1] = sub->vdim[1];
      }
      cov->deterministic &= sub->deterministic;
    };
    //  printf("C OK err=%d\n", err);

    if (i==0) {
      cov->vdim[0]=sub->vdim[0];  // to do: inkonsistent mit vorigen Zeilen !!
      cov->vdim[1]=sub->vdim[1];  // to do: inkonsistent mit vorigen Zeilen !!
      if (cov->vdim[0] <= 0) BUG;
      cov->matrix_indep_of_x = sub->matrix_indep_of_x;
    } else {
      cov->matrix_indep_of_x &= sub->matrix_indep_of_x;
      if (cov->vdim[0] != sub->vdim[0] || cov->vdim[1] != sub->vdim[1]) {
	SERR4("multivariate dimensionality is different in the submodels (%s is %d-variate; %s is %d-variate)", NICK(cov->sub[0]), cov->vdim[0], NICK(sub), sub->vdim[0]);
      }
    }
  }

  // !! incorrect  !!
  // cov->semiseparatelast = false; 
  //cov->separatelast = false; 

  //  printf("OK err=%d ende\n", err);
   //  PMI(cov, -1);

  return NOERROR;
}





// see private/old.code/operator.cc for some code including variable locations
void select(double *x, cov_model *cov, double *v) {
  int len,
    *element = PINT(SELECT_SUBNR);
  cov_model *sub = cov->sub[*element];
  assert(cov->vdim[0] == cov->vdim[1]);
  if (*element >= cov->nsub) error("select: element out of range");
  COV(x, sub, v);
  if ( (len = cov->nrow[SELECT_SUBNR]) > 1) {
    int i, m,
      vdim = cov->vdim[0],
      vsq = vdim * vdim;
    ALLOC_EXTRA(z, vsq);
    for (i=1; i<len; i++) {
      sub = cov->sub[element[i]];
      COV(x, sub, z);
      for (m=0; m<vsq; m++) v[m] += z[m]; 
    }
  }
}
  

void covmatrix_select(cov_model *cov, double *v, int *nonzeros) {
  int len = cov->nrow[SELECT_SUBNR];
  
  if (len == 1) {
    int element = P0INT(SELECT_SUBNR);
    cov_model *next = cov->sub[element];
    if (element >= cov->nsub) error("select: element out of range");
    CovList[next->nr].covmatrix(next, v, nonzeros);

    // {  int i,j,k, tot=Loc(cov)->totalpoints; printf("\nXcovmat select\n");
    //  for (k=i=0; i<tot*tot; i+=tot) {
    //   for (j=0; j<tot; j++) printf("%f ", v[k++]);
    //  printf("\n");  }}

    // crash();
    //  PMI(next);
      
  }  else StandardCovMatrix(cov, v, nonzeros);
}

char iscovmatrix_select(cov_model VARIABLE_IS_NOT_USED *cov) {  return 2; }

int checkselect(cov_model *cov) {
  int err;

  assert(cov->Splus == NULL);

  if (!isCartesian(cov->isoown)) return ERRORNOTCARTESIAN;
  kdefault(cov, SELECT_SUBNR, 0);
  if ((err = checkplus(cov)) != NOERROR) return err;

  if ((err = checkkappas(cov)) != NOERROR) return err;

  EXTRA_STORAGE;

  return NOERROR;
}


void rangeselect(cov_model VARIABLE_IS_NOT_USED *cov, range_type *range){
  range->min[SELECT_SUBNR] = 0;
  range->max[SELECT_SUBNR] = MAXSUB-1;
  range->pmin[SELECT_SUBNR] = 0;
  range->pmax[SELECT_SUBNR] = MAXSUB-1;
  range->openmin[SELECT_SUBNR] = false;
  range->openmax[SELECT_SUBNR] = false;
}




void plusStat(double *x, cov_model *cov, double *v){
  cov_model *sub;
  int i, m,
    nsub=cov->nsub,
    vdim = cov->vdim[0],
    vsq = vdim * vdim;
  assert(cov->vdim[0] == cov->vdim[1]);
  ALLOC_EXTRA(z, vsq);
  
  for (m=0; m<vsq; m++) v[m] = 0.0;
  for(i=0; i<nsub; i++) {
    sub = cov->sub[i];
    if (cov->typus == sub->typus) {
      COV(x, sub, z);
      for (m=0; m<vsq; m++) v[m] += z[m]; 
    }
  }
}

void plusNonStat(double *x, double *y, cov_model *cov, double *v){
  cov_model *sub;
  int i, m,
    nsub=cov->nsub,
    vdim = cov->vdim[0],
    vsq = vdim * vdim;
  assert(cov->vdim[0] == cov->vdim[1]);
  ALLOC_EXTRA(z, vsq);
  for (m=0; m<vsq; m++) v[m] = 0.0;
  for(i=0; i<nsub; i++) {
    sub = cov->sub[i];
    if (cov->typus == sub->typus) {
      NONSTATCOV(x, y, sub, z);
      for (m=0; m<vsq; m++) v[m] += z[m]; 
    }

    //printf("i=%d %f %f\n", i, v[0], z[0]);
  }
  // printf("plus nonstat x=%f %f\n", *x, *v);
}

void Dplus(double *x, cov_model *cov, double *v){
  cov_model *sub;
  double v1;
  int n = cov->nsub, i;
  *v = 0.0;
  for (i=0; i<n; i++) { 
    sub = cov->sub[i];
    if (cov->typus == sub->typus) {
      Abl1(x, sub, &v1);
      (*v) += v1;
    }
  }
}

void DDplus(double *x, cov_model *cov, double *v){
  cov_model *sub;
  double v1;
  int n = cov->nsub, i;
  *v = 0.0;
  for (i=0; i<n; i++) { 
    sub = cov->sub[i];
    if (cov->typus == sub->typus) {
      Abl2(x, sub, &v1);
      (*v) += v1;
    }
  }
}


int checkplus(cov_model *cov) {
  int err, i;
  if ((err = checkplusmal(cov)) != NOERROR) {
    return err;
  }
  
  if (cov->domown == DOMAIN_MISMATCH) return ERRORNOVARIOGRAM;
  if (cov->nsub == 0) cov->pref[SpectralTBM] = PREF_NONE;

  if (isPosDef(cov) && cov->domown == XONLY) cov->logspeed = 0.0;
  else if (isNegDef(cov) && cov->domown == XONLY) {
    cov->logspeed = 0.0;
    for (i=0; i<cov->nsub; i++) {
      cov_model *sub = cov->sub[i];
      if (cov->typus == sub->typus) {
	if (ISNAN(sub->logspeed)) {
	  cov->logspeed = RF_NA;
	  break;
	} else cov->logspeed += sub->logspeed;
      }
    }
  } else cov->logspeed = RF_NA;

  EXTRA_STORAGE;

  return NOERROR;

  // spectral mit "+" funktioniert, falls alle varianzen gleich sind,
  // d.h. nachfolgend DOLLAR die Varianzen gleich sind oder DOLLAR nicht
  // auftritt; dann zufaellige Auswahl unter "+"
}


bool Typeplus(Types required, cov_model *cov, int depth) {
  // assert(false);
  bool allowed = TypeConsistency(ShapeType, required) || required==TrendType;
    //||required==ProcessType ||required==GaussMethodType; not yet allowed;to do

  // printf("allowed %d\n", allowed);
  if (!allowed) return false;
  int i;
  for (i=0; i<cov->nsub; i++) {
    //printf("type plus %s %s\n",TYPENAMES[required], NICK(cov->sub[i]));
    if (TypeConsistency(required, cov->sub[i], depth-1)) return true;
  }  
  // printf("none\n");
  
  return false;
}

void spectralplus(cov_model *cov, gen_storage *s, double *e){
  int nr;
  double dummy;
  spec_properties *cs = &(s->spec);
  double *sd_cum = cs->sub_sd_cum;

  nr = cov->nsub - 1;
  dummy = UNIFORM_RANDOM * sd_cum[nr];
  if (ISNAN(dummy)) BUG;
  while (nr > 0 && dummy <= sd_cum[nr - 1]) nr--;
  cov_model *sub = cov->sub[nr];
  SPECTRAL(sub, s, e);  // nicht gatternr
}


int structplus(cov_model *cov, cov_model VARIABLE_IS_NOT_USED **newmodel){
  int m, err;  
  switch(cov->role) {
  case ROLE_COV :  return NOERROR;
  case ROLE_GAUSS : 
    if (isProcess(cov->typus)) {
      //assert(cov->calling != NULL && (isInterface(cov->calling->typus) ||
      ///				      isProcess(cov->calling->typus)));
      assert(cov->nr != PLUS_PROC);
      assert(cov->nr == PLUS);
      //cov->nr = PLUS_PROC;
      BUG;
      //return structplusproc(cov, newmodel); // kein S-TRUCT !!
    }
    if (cov->Splus != NULL) BUG;
    for (m=0; m<cov->nsub; m++) {
      cov_model *sub = cov->sub[m];
      if ((err = STRUCT(sub, newmodel))  > NOERROR) {
      //	PMI(cov->Splus->keys[m]);
	//	assert(false);
	//printf("end plus\n");
	return err;
      }
    }
    break;
  default :
    SERR2("role '%s' not allowed for '%s'", ROLENAMES[cov->role],
	  NICK(cov));
  }
  return NOERROR;
}


int initplus(cov_model *cov, gen_storage *s){
  int i, err,
    vdim = cov->vdim[0];
  if (cov->vdim[0] != cov->vdim[1]) BUG;

  for (i=0; i<vdim; i++) cov->mpp.maxheights[i] = RF_NA;

  if (cov->role == ROLE_GAUSS) {
    spec_properties *cs = &(s->spec);
    double *sd_cum = cs->sub_sd_cum;
 
    for (i=0; i<cov->nsub; i++) {
      cov_model *sub = cov->Splus == NULL ? cov->sub[i] : cov->Splus->keys[i];

      //print("initplus %d %d initialised=%d\n",i,cov->nsub, sub->initialised);
      //PMI(sub);
      
      if (sub->pref[Nothing] > PREF_NONE) { // to do ??
	// for spectral plus only
	COV(ZERO, sub, sd_cum + i);
	if (i>0) sd_cum[i] += sd_cum[i-1];
      }
      cov->sub[i]->Sgen = (gen_storage *) MALLOC(sizeof(gen_storage));
      if ((err = INIT(sub, cov->mpp.moments, s)) != NOERROR) {
	//  AERR(err);
	return err;
      }
      sub->simu.active = true;
    }

    // assert(false);
  
    cov->fieldreturn = cov->Splus != NULL;
    cov->origrf = false;
    if (cov->Splus != NULL) cov->rf = cov->Splus->keys[0]->rf;
     
    return NOERROR;
  }

  /*
    pref_type pref;

    for (m=0; m<Forbidden; m++) pref[m] = PREF_BEST;
    if (meth->cov->user[0] == 0 || meth->cov->user[1] == 0) {
 
    // i.e. user defined
    pref_type pref;
    pref[CircEmbed] = pref[TBM] = pref[Direct] = pref[Sequential] 
    = PREF_NONE;
    for (m=0; m<Nothing; m++) { // korrekt auch fuer MaxStable?
    //	  print("%d %d\n", m, pref[m]);
    //	  print("%d %d\n", meth->cov->user[m]);
    if (pref[m] > 0 &&  meth->cov->user[m] > 0) {
    break;
    }
    }
    if (m == Nothing) return ERRORSUBMETHODFAILED;
    }
  */
     
  else if (cov->role == ROLE_COV) {    
    return NOERROR;
  }

  return ERRORFAILED;
}


void doplus(cov_model *cov, gen_storage *s) {
  int i;
  
  if (cov->role == ROLE_GAUSS && cov->method==SpectralTBM) {
    ERR("error in doplus with spectral");
  }
  
  for (i=0; i<cov->nsub; i++) {
    cov_model *sub = cov->Splus==NULL ? cov->sub[i] : cov->Splus->keys[i];
    DO(sub, s);
  }
}




void covmatrix_plus(cov_model *cov, double *v, int *nonzeros) {
  location_type *loc = Loc(cov);
  //  cov_fct *C = CovList + cov->nr; // nicht gatternr
  long i, 
    totalpoints = loc->totalpoints,
    vdimtot = totalpoints * cov->vdim[0],
    vdimtotSq = vdimtot * vdimtot;
  int
    nsub = cov->nsub;
  bool is = iscovmatrix_plus(cov) >= 2;
  double *mem = NULL;
  if (is && nsub > 1) {
    ALLOC_EXTRA2(MEM, vdimtotSq);
    mem = MEM;
    is = mem != NULL;
  }
  
  if (is) {
    // APMI(cov);

    //cov_model *sub = cov->sub[0];
    int j;
    if (PisNULL(SELECT_SUBNR)) PALLOC(SELECT_SUBNR, 1, 1);
    P(SELECT_SUBNR)[0] = 0;
    CovList[SELECT].covmatrix(cov, v, nonzeros);
    for (i=1; i<nsub; i++) {
      if (Loc(cov->sub[i])->totalpoints != totalpoints) BUG;
      P(SELECT_SUBNR)[0] = i;
      CovList[SELECT].covmatrix(cov, mem, nonzeros);
      for (j=0; j<vdimtotSq; j++) v[j] += mem[j];
    }
    int nz = 0;
    for (j=0; j<vdimtotSq; nz += v[j++] != 0.0);
    *nonzeros = nz;
  } else StandardCovMatrix(cov, v, nonzeros);  
}

char iscovmatrix_plus(cov_model *cov) {
  char max=0, is;
  int i,
    nsub = cov->nsub;
  for (i=0; i<nsub; i++) {
    cov_model *sub = cov->sub[i];
    is = CovList[sub->nr].is_covmatrix(sub);
    if (is > max) max = is;
  }
  return max;
}


void malStat(double *x, cov_model *cov, double *v){
  cov_model *sub;
  int i, m,
    nsub=cov->nsub,
    vdim = cov->vdim[0],
    vsq = vdim * vdim;
  ALLOC_EXTRA(z, vsq);
  
  assert(x[0] >= 0.0 || cov->xdimown > 1);
  for (m=0; m<vsq; m++) v[m] = 1.0;
  for(i=0; i<nsub; i++) {
    sub = cov->sub[i];
    COV(x, sub, z);
    for (m=0; m<vsq; m++) v[m] *= z[m]; 
  }
}

void logmalStat(double *x, cov_model *cov, double *v, double *sign){
  cov_model *sub;
  int i, m,
    nsub=cov->nsub,
    vdim = cov->vdim[0],
    vsq = vdim * vdim;
  ALLOC_EXTRA(z, vsq);
  ALLOC_EXTRA(zsign, vsq);

  assert(cov->vdim[0] == cov->vdim[1]); 
  assert(x[0] >= 0.0 || cov->xdimown > 1);
  for (m=0; m<vsq; m++) {v[m] = 0.0; sign[m]=1.0;}
  for(i=0; i<nsub; i++) {
    sub = cov->sub[i];
    LOGCOV(x, sub, z, zsign);
    for (m=0; m<vsq; m++) {
      v[m] += z[m]; 
      sign[m] *= zsign[m];
    }
  }
}

void malNonStat(double *x, double *y, cov_model *cov, double *v){
  cov_model *sub;
  int i, m, nsub=cov->nsub,
    vdim = cov->vdim[0],
    vsq = vdim * vdim;
  ALLOC_EXTRA(z, vsq);

  assert(cov->vdim[0] == cov->vdim[1]);
  assert(cov->vdim[0] == cov->vdim[1]);

  for (m=0; m<vsq; m++) v[m] = 1.0;
  for(i=0; i<nsub; i++) {
    sub = cov->sub[i];
    NONSTATCOV(x, y, sub, z);
    for (m=0; m<vsq; m++) v[m] *= z[m]; 
  }
}

void logmalNonStat(double *x, double *y, cov_model *cov, double *v, 
		   double *sign){
  cov_model *sub;
  int i, m, nsub=cov->nsub,
    vdim = cov->vdim[0],
    vsq = vdim * vdim;
  assert(cov->vdim[0] == cov->vdim[1]);
  ALLOC_EXTRA(z, vsq);
  ALLOC_EXTRA(zsign, vsq);
  for (m=0; m<vsq; m++) {v[m] = 0.0; sign[m]=1.0;}
  for(i=0; i<nsub; i++) {
    sub = cov->sub[i];
    LOGNONSTATCOV(x, y, sub, z, zsign);
    for (m=0; m<vsq; m++) {
      v[m] += z[m]; 
      sign[m] *= zsign[m];
    }
  }
}

void Dmal(double *x, cov_model *cov, double *v){
  cov_model *sub;
  double c[MAXSUB], d[MAXSUB];
  int n = cov->nsub, i;
  for (i=0; i<n; i++) {
    sub = cov->sub[i];
    COV(x, sub, c + i);
    Abl1(x, sub, d + i);
  }
  *v = 0.0;
  for (i=0; i<n; i++) {
    double zw = d[i];
    int j;
    for (j=0; j<n; j++) if (j!=i) zw *= c[j]; 
    (*v) += zw; 
  }
}
int checkmal(cov_model *cov) {
  cov_model *next1 = cov->sub[0];
  cov_model *next2 = cov->sub[1];
  int err;

  if (next2 == NULL) next2 = next1;
  
  if ((err = checkplusmal(cov)) != NOERROR) return err;

  if (cov->domown == DOMAIN_MISMATCH || !isPosDef(cov)) {
    return ERRORNOVARIOGRAM;
  }
  // to do istcftype und allgemeiner typ zulassen

  cov->logspeed = cov->domown == XONLY ? 0 : RF_NA;
  
  if (cov->xdimown >= 2) cov->pref[TBM] = PREF_NONE;
  if (cov->xdimown==2 && cov->nsub == 2 && 
      isAnyDollar(next1) && isAnyDollar(next2)) {
    double *aniso1 = PARAM(next1, DANISO),
      *aniso2 = PARAM(next2, DANISO);
    if (aniso1 != NULL && aniso2 != NULL) {
      if (aniso1[0] == 0.0 && next1->ncol[DANISO] == 1) {
	cov->pref[TBM] = next2->pref[TBM];
      } else if (aniso2[0] == 0.0 && next2->ncol[DANISO] == 1) {
	cov->pref[TBM] = next1->pref[TBM];
      }
    }
  }
  
  EXTRA_STORAGE;

  return NOERROR;
}

bool Typemal(Types required, cov_model *cov, int depth) {
  if (!isShape(required)) return false;
  int i;
  for (i=0; i<cov->nsub; i++) {
    if (!TypeConsistency(required, cov->sub[i], depth-1)) return false;
  }
  return true;
}


int initmal(cov_model *cov, gen_storage VARIABLE_IS_NOT_USED *s){
//  int err;
//  return err;
  return ERRORFAILED;
  int i, 
    vdim = cov->vdim[0];
  if (cov->vdim[0] != cov->vdim[1]) BUG;

  for (i=0; i<vdim; i++) cov->mpp.maxheights[i] = RF_NA;
}
void domal(cov_model VARIABLE_IS_NOT_USED *cov, gen_storage VARIABLE_IS_NOT_USED *s){
  BUG;
}


//////////////////////////////////////////////////////////////////////
// mpp-plus


#define PLUS_P 0 // parameter
int  CheckAndSetP(cov_model *cov){   
  int n,
    nsub = cov->nsub;
  double 
    cump = 0.0;
  if (PisNULL(PLUS_P)) {
    assert(cov->nrow[PLUS_P] == 0 && cov->ncol[PLUS_P] == 0);
    PALLOC(PLUS_P, nsub, 1);
    for (n=0; n<nsub; n++) P(PLUS_P)[n] = 1.0 / (double) nsub;    
  } else {
    cump = 0.0;
    for(n = 0; n<nsub; n++) {
      cump += P(PLUS_P)[n];
      //printf("cump %f %f\n",  cump , P(PLUS_P)[n]);
      if (cump > 1.0 && n+1<nsub) return ERRORATOMP; 
    }
    if (cump != 1.0) {
      if (nsub == 1) {
	warning("the p-values do not sum up to 1.\nHere only one p-value is given which must be 1.0");
	P(PLUS_P)[0] = 1.0;
      } else {
	//printf("%e\n", 1-P(PLUS_P)[nsub-1] );
	if (cump < 1.0 && P(PLUS_P)[nsub-1] == 0) {
	  WARNING1("The value of the last component of '%s' is increased.",
		   KNAME(PLUS_P)); 

	}
	else SERR1("The components of '%s' do not sum up to 1.", KNAME(PLUS_P));
	P(PLUS_P)[nsub-1] = 1.0 - (cump - P(PLUS_P)[nsub-1]);
      }  
    }
  }
  return NOERROR;
}

void kappamppplus(int i, cov_model *cov, int *nr, int *nc){
  *nr = cov->nsub;
  *nc = i < CovList[cov->nr].kappas ? 1 : -1;
}

void mppplus(double *x, cov_model *cov, double *v) { 
  int i, n,
    vdim = cov->vdim[0],
    vdimSq = vdim * vdim;
  cov_model *sub;
  ALLOC_EXTRA(z, vdimSq);

  if (cov->role == ROLE_COV) {  
    for (i=0; i<vdimSq; i++) v[i] = 0.0;
    for (n=0; n<cov->nsub; n++, sub++) {
      sub = cov->sub[n];
      COV(x, sub, z); // urspruenglich : covlist[sub].cov(x, cov, z); ?!
      for (i=0; i<vdimSq; i++) v[i] += P(PLUS_P)[n] * z[i];
    }
  } else {
    assert(hasPoissonRole(cov));
    BUG;
  }
}

int checkmppplus(cov_model *cov) { 
  int err, 
    size = 1;
  cov->maxdim = MAXMPPDIM;
  
  if ((err = checkplusmal(cov)) != NOERROR) {
    // printf("checkmppplus error %d %s\n", err, ERRORSTRING);
    return err;
  }

  if ((err = CheckAndSetP(cov)) != NOERROR) return(err);

  if (cov->q == NULL) QALLOC(size);
  EXTRA_STORAGE;
  
  return NOERROR;
}



void rangempplus(cov_model VARIABLE_IS_NOT_USED *cov, range_type *range){ 
  range->min[PLUS_P] = 0.0;
  range->max[PLUS_P] = 1.0;
  range->pmin[PLUS_P] = 0.0;
  range->pmax[PLUS_P] = 1.0;
  range->openmin[PLUS_P] = false;
  range->openmax[PLUS_P] = false;
}

int struct_mppplus(cov_model *cov, cov_model **newmodel){
  int m,
    //nsub = cov->nsub,
    err = NOERROR;

  //if (cov->calling == NULL || cov->calling->ownloc == NULL) 
  // SERR("mppplus does not seem to be used in the right context");
  
  if (!hasMaxStableRole(cov) && !hasPoissonRole(cov)) {
    SERR("method is not based on Poisson point process");
  }


  return ERRORNOTPROGRAMMEDYET;

  
  // Ausnahme: mppplus wird separat behandelt:
  // if (nr == MPPPLUS) return S TRUCT(shape, NULL);
 
  ASSERT_NEWMODEL_NOT_NULL;
  NEW_STORAGE(plus);
  plus_storage *s = cov->Splus;

  for (m=0; m<cov->nsub; m++) {
    cov_model *sub = cov->sub[m];          
    if (s->keys[m] != NULL) COV_DELETE(s->keys + m);    
    if ((err = covcpy(s->keys + m, sub)) != NOERROR) return err;
    if ((err = addShapeFct(s->keys + m)) != NOERROR) return err;
    s->keys[m]->calling = cov;
  }
  return NOERROR;
}


int init_mppplus(cov_model *cov, gen_storage *S) {
  cov_model  *sub;
  double M2[MAXMPPVDIM], M2plus[MAXMPPVDIM], Eplus[MAXMPPVDIM], 
    maxheight[MAXMPPVDIM];
  int i,n, err,
    vdim = cov->vdim[0];
  if (cov->vdim[0] != cov->vdim[1]) BUG;
  if (vdim > MAXMPPVDIM) BUG;
  ext_bool 
    loggiven = SUBMODEL_DEP,
    fieldreturn = SUBMODEL_DEP;
  pgs_storage *pgs = NULL;
  for (i=0; i<vdim; i++) {
    maxheight[i] = RF_NEGINF;
    M2[i] = M2plus[i] = Eplus[i] = 0.0;
  }
    
  NEW_STORAGE(pgs);
  pgs = cov->Spgs;
  pgs->totalmass = 0.0;
  
  for (n=0; n<cov->nsub; n++) {
    sub = cov->sub[n];
    //if (!sub->mpp.loc_done) 
    //SERR1("submodel %d of '++': the location of the shapes is not defined", 
    // n);
    if ((err = INIT(sub, cov->mpp.moments, S)) != NOERROR) return err;
    //if (!sub->mpp.loc_done) SERR("locations are not initialised");
    if (n==0) {
      loggiven = sub->loggiven;
      fieldreturn = sub->fieldreturn;
    } else {
      if (loggiven != sub->loggiven) loggiven = SUBMODEL_DEP;
      if (fieldreturn != sub->loggiven) fieldreturn = SUBMODEL_DEP;
    }
    pgs->totalmass += sub->Spgs->totalmass * P(PLUS_P)[n];
    for (i=0; i<vdim; i++)   
      if (cov->mpp.maxheights[i] > maxheight[i]) 
	maxheight[i] = cov->mpp.maxheights[i];
    loggiven &= cov->loggiven;

    // Achtung cov->mpp.mM2 und Eplus koennten nicht direkt gesetzt
    // werden, da sie vom CovList[sub->nr].Init ueberschrieben werden !!
    
    if (cov->mpp.moments >= 1) {
      int nmP1 = sub->mpp.moments + 1;
      for (i=0; i<vdim; i++) {
	int idx = i * nmP1;
	Eplus[i] += PARAM0(sub, PLUS_P) * sub->mpp.mMplus[idx + 1]; 
      }
      if (cov->mpp.moments >= 2) {
	for (i=0; i<vdim; i++) {
	  int idx = i * nmP1;
	  M2[i] += PARAM0(sub, PLUS_P)  * sub->mpp.mM[idx + 2];
	  M2plus[i] += PARAM0(sub, PLUS_P) * sub->mpp.mM[idx + 2];
	}
      }
    }
    //assert(sub->mpp.loc_done);
  }
  for (i=0; i<vdim; i++) cov->mpp.maxheights[i] = maxheight[i];
  //cov->mpp.refsd = RF_NA;
  //  cov->mpp.refradius = RF_NA;

  
  if (cov->mpp.moments >= 1) {
    int nmP1 = cov->mpp.moments + 1;
    for (i=0; i<vdim; i++) {
      int idx = i * nmP1;
      cov->mpp.mMplus[idx + 1] = Eplus[i];
      cov->mpp.mM[idx + 1] = RF_NA;
    }
    if (cov->mpp.moments >= 2) {
      for (i=0; i<vdim; i++) {
	 int idx = i * nmP1;
	 cov->mpp.mM[idx + 2] = M2[i];
	 cov->mpp.mMplus[idx + 2] = M2plus[i];
      }
    }
  }
  
  cov->loggiven = loggiven;
  cov->fieldreturn = fieldreturn;
  //cov->mpp.loc_done = true;       
  cov->origrf = false;
  cov->rf = NULL;

  return NOERROR;
}

void do_mppplus(cov_model *cov, gen_storage *s) {
  cov_model *sub;
  double subselect = UNIFORM_RANDOM;
  int i, subnr,
    vdim = cov->vdim[0];
  assert(cov->vdim[0] == cov->vdim[1]);
  for (subnr=0; (subselect -= PARAM0(cov->sub[subnr], PLUS_P)) > 0; subnr++);
  cov->q[0] = (double) subnr;
  sub = cov->sub[subnr];
  
  CovList[sub->nr].Do(sub, s);  // nicht gatternr
  for (i=0; i<vdim; i++)   
    cov->mpp.maxheights[i] = sub->mpp.maxheights[i];
  cov->fieldreturn = sub->fieldreturn;
  cov->loggiven = sub->loggiven;
  // PrintMPPInfo(cov, s);
}

//////////////////////////////////////////////////////////////////////
// PROCESSES
//////////////////////////////////////////////////////////////////////


int structSproc(cov_model *cov, cov_model **newmodel) {
  cov_model
    *next = cov->sub[DOLLAR_SUB],
    *Aniso = cov->kappasub[DAUSER];
  int dim, err; 
  // cov_model *sub;

  //   printf("%d\n", cov->role); assert(cov->role == ROLE_GAUSS);
  
  //   assert(false);

  //    assert(false);
 
   assert(isDollarProc(cov));
   //printf("her %d %d\n", cov->role, ROLE_GAUSS);

  if (Aniso != NULL && !Aniso->deterministic) {
     SERR1("complicated models including arbitrary functions for '%s' cannot be simulated yet", KNAME(DAUSER));
  }



  switch (cov->role) {
  case ROLE_GAUSS :
    ASSERT_NEWMODEL_NULL;
    if (cov->key != NULL) COV_DELETE(&(cov->key));

    if (cov->prevloc->distances) 
      SERR("distances do not allow for more sophisticated simulation methods");
    
    if (Aniso!= NULL) {
      //A
      // crash();
      Transform2NoGrid(cov, false, true); 
      
      location_type *loc = Loc(cov);
      dim = loc->timespacedim;
      int 
	bytes = dim * sizeof(double);
      long i,
	total = loc->totalpoints;
      if (dim != Aniso->vdim[0]) BUG;
      double *v = NULL,
	*x = loc->x;
      assert(x != NULL);
      assert(!loc->grid);
      if ((v = (double*) MALLOC(bytes)) == NULL) return ERRORMEMORYALLOCATION;
      for (i=0; i<total; i++, x+=dim) {
	FCTN(x, Aniso, v);
	//printf("x = %f %f  v=%f %f %d\n", x[0], x[1], v[0], v[1], bytes);
	MEMCOPY(x, v, bytes);
      }
      free(v);
      v = NULL;
    } else {
      Transform2NoGrid(cov, true, false); 
    }

    //assert(false);
    // printf("%f %f %f\n",
	   //	   Loc(cov)->xgr[0][0], Loc(cov)->xgr[0][1], Loc(cov)->xgr[0][2]);
    // assert(false);
 
    if ((err = covcpy(&(cov->key), next)) != NOERROR) return err;

    if (!isGaussProcess(next)) addModel(&(cov->key), GAUSSPROC);
    SetLoc2NewLoc(cov->key, Loc(cov));
    //   APMI(cov);assert(false);
    
    cov_model *key;
    key = cov->key;
    assert(key->calling == cov);    
    
    dim = key->prevloc->timespacedim;
    if ((err = CHECK(key, dim, dim, ProcessType, XONLY, CARTESIAN_COORD,
		     cov->vdim, cov->role)) != NOERROR) return err;
    err = STRUCT(cov->key, NULL);
    //    cov->initialised = err==NOERROR && key->initialised;

    //    APMI(cov);

    return err;
  default :
    //  PMI(cov, "structS");
    SERR2("%s: changes in scale/variance not programmed yet for '%s'", 
	  NICK(cov), ROLENAMES[cov->role]);      
  }
     
  return NOERROR;
}




int initSproc(cov_model *cov, gen_storage *s){
  // am liebsten wuerde ich hier die Koordinaten transformieren;
  // zu grosser Nachteil ist dass GetDiameter nach trafo 
  // grid def nicht mehr ausnutzen kann -- umgehbar?!

  //cov_model *next = cov->sub[DOLLAR_SUB];
  //mppinfotype *info = &(s->mppinfo);
  //  location_type *loc = cov->prevloc;
  int 
    err = NOERROR;

  //  assert(false);
  
  cov_model 
    *key = cov->key;
    //*sub = key == NULL ? next : key;
  location_type *prevloc = cov->prevloc;

  assert(key != NULL);
  
  if ((err = INIT(key, 0, s)) != NOERROR) {
    return err;
  }
  
  key->simu.active = true; 
  assert(s != NULL);

  cov->fieldreturn = true;

  if ((cov->origrf = cov->ownloc != NULL &&
       cov->ownloc->totalpoints != prevloc->totalpoints)) {
    assert(prevloc->grid);
    int dim = prevloc->timespacedim;
    if (cov->vdim[0] != cov->vdim[1]) BUG;
    cov->rf = (res_type*) MALLOC(sizeof(res_type) *
				 cov->vdim[0] * 
				 prevloc->totalpoints);
    DOLLAR_STORAGE;
 
    int d,
      *proj = PINT(DPROJ),
      bytes = dim * sizeof(int),
      *cumsum = cov->Sdollar->cumsum = (int*) MALLOC(bytes),
      *total = cov->Sdollar->total = (int*) MALLOC(bytes),
      *len = cov->Sdollar->len = (int*) MALLOC(bytes);     
    cov->Sdollar->nx = (int*) MALLOC(bytes); 
    
    for (d=0; d<dim; d++) {
      cumsum[d] = 0;
      len[d] = prevloc->xgr[d][XLENGTH];
    }
    if (proj != NULL) {
      int 
	nproj = cov->nrow[DPROJ];
      d = 0;
      cumsum[proj[d] - 1] = 1;
      for (d = 1; d < nproj; d++) {
	cumsum[proj[d] - 1] =
	  cumsum[proj[d - 1] - 1] * prevloc->xgr[d - 1][XLENGTH];
      }
      for (d=0; d<dim;d++) 
	total[d] = cumsum[d] * prevloc->xgr[d][XLENGTH];
    } else {
      int i,
	iold = 0,
	nrow = cov->nrow[DANISO],
	ncol = cov->ncol[DANISO];
      double *A = P(DANISO);
      for (d=0; d<ncol; d++, A += nrow) {
	for (i = 0; i < nrow && A[i] == 0.0; i++);
	if (i == nrow) i = nrow - 1;
	if (d > 0) {
	  cumsum[i] = cumsum[iold] * prevloc->xgr[d - 1][XLENGTH];
	} else { // d ==0
	  cumsum[i] = 1;
	}
	iold = i;
	for (i++; i < nrow; i++) if (A[i] != 0.0) BUG;  // just a check
      }
    }
  } else {
    cov->rf = cov->key->rf;
  }
 
  //  PMI(cov,-1); assert(false);
 
  return NOERROR;
}


void doSproc(cov_model *cov, gen_storage *s){

  if (hasMaxStableRole(cov) || hasPoissonRole(cov)) {
    cov_model *next = cov->sub[DOLLAR_SUB];
    
    cov_model
      *varM = cov->kappasub[DVAR],
      *scaleM = cov->kappasub[DSCALE];
 
    int i,
      vdim = cov->vdim[0];
    assert(cov->vdim[0] == cov->vdim[1]);

    if (varM != NULL && !varM->deterministic) {
      assert(!PisNULL(DVAR));
      VTLG_R(NULL, varM, P(DVAR));
    }

    if (scaleM != NULL && !scaleM->deterministic) {
      assert(!PisNULL(DSCALE));
      VTLG_R(NULL, scaleM, P(DSCALE));
    }
    
    DO(next, s);// nicht gatternr
    for (i=0; i<vdim; i++)   
      cov->mpp.maxheights[i] = next->mpp.maxheights[i] * P0(DVAR);

  }

  else if (cov->role == ROLE_GAUSS) {    
    assert(cov->key != NULL);
    double 
      *res = cov->key->rf,
      sd = sqrt(P0(DVAR));
    int i,
      totalpoints = Gettotalpoints(cov);

    DO(cov->key, s); 

    if (sd != 1.0) 
      for (i=0; i<totalpoints; i++) {
	res[i] *= (res_type) sd;
    }

  }
  
  else BUG;
 
  if (cov->origrf) {
    assert(cov->prevloc->grid);
    int i, zaehler, d,
      dim = cov->prevloc->timespacedim,
      *cumsum = cov->Sdollar->cumsum,
      *nx = cov->Sdollar->nx,
      *len = cov->Sdollar->len,
      *total = cov->Sdollar->total;
    assert(cov->key != NULL);
    res_type
      *res = cov->rf,
      *rf = cov->key->rf;
    
    assert(nx != NULL && total != NULL && cumsum != NULL);
    for (d=0; d<dim; d++) {
      nx[d] = 0;
    }
    zaehler = 0;
    i = 0;      
    
    while (true) {
      res[i++] = rf[zaehler];
      // printf("%f %d %d", res[i], zaehler,i);
      d = 0;			
      nx[d]++;			
      zaehler += cumsum[d];	
      while (nx[d] >= len[d]) {	
	nx[d] = 0;		
	zaehler -= total[d];
	if (++d >= dim) break;	
	nx[d]++;			
	zaehler += cumsum[d];					
      }
      if (d >= dim) break;			
    }
  }

}


int checkplusmalproc(cov_model *cov) {
  cov_model *sub;
  int i, err, 
    dim = cov->tsdim, 
    xdim = cov->xdimown,
    role = cov->role;
  Types type = ProcessType; 
  domain_type dom = cov->domown;
  isotropy_type iso = cov->isoown;

  //PMI(cov);
  assert(cov->Splus != NULL);

  for (i=0; i<cov->nsub; i++) {

    sub = cov->Splus->keys[i];
     
    if (sub == NULL) 
      SERR("named submodels are not given in a sequence.");

    if (!TypeConsistency(type, sub, 0)) return ERRORTYPECONSISTENCY;
    if ((err= CHECK(sub, dim, xdim, type, dom, iso, SUBMODEL_DEP, role))
	!= NOERROR) return err;

   if (i==0) {
      cov->vdim[0]=sub->vdim[0];  // to do: inkonsistent mit vorigen Zeilen !!
      cov->vdim[1]=sub->vdim[1];  // to do: inkonsistent mit vorigen Zeilen !!
    } else {
      if (cov->vdim[0] != sub->vdim[0] || cov->vdim[1] != sub->vdim[1]) {
	SERR("multivariate dimensionality must be equal in the submodels");
      }
    }

    //    printf("OK plusprocess\n");
  }

  return NOERROR;
}



int checkplusproc(cov_model *cov) {
  int err;
  if ((err = checkplusmalproc(cov)) != NOERROR) {
    return err;
  }
  return NOERROR;
}


int structplusmalproc(cov_model *cov, cov_model VARIABLE_IS_NOT_USED**newmodel){
  int m, err;

  switch(cov->role) {
  case ROLE_GAUSS : 
    {
      location_type *loc = Loc(cov);
      NEW_STORAGE(plus);
      plus_storage *s =cov->Splus;
      for (m=0; m<cov->nsub; m++) {

	cov_model *sub = cov->sub[m];

	if (s->keys[m] != NULL) COV_DELETE(s->keys + m);
	if ((err =  covcpy(s->keys + m, sub)) != NOERROR) {
	  return err;
	}
	assert(s->keys[m] != NULL);
	assert(s->keys[m]->calling == cov);
	
	if (PL >= PL_STRUCTURE) {
	  LPRINT("plus: trying initialisation of submodel #%d (%s).\n", m+1, 
		 NICK(sub));
	}
	
	addModel(s->keys + m, GAUSSPROC);
	s->keys[m]->calling = cov;
	//	cov_model *fst = cov; while (fst->calling != NULL) fst = fst->calling; 

	//	assert(false);

	err = CHECK(s->keys[m], loc->timespacedim, loc->timespacedim,
		    ProcessType, XONLY, CARTESIAN_COORD, cov->vdim, ROLE_GAUSS);
	if (err != NOERROR) {
	  //	  
	  return err;
	}
	//APMI(s->keys[m]);
	
	//if (m==1) APMI(s->keys[m]);
	if ((s->struct_err[m] =
	     err = STRUCT(s->keys[m], NULL))  > NOERROR) {
	  //	PMI(s->keys[m]);
	  //	assert(false);
	  //printf("end plus\n");
	  //	  	  AERR(err);
	  return err;
	}
	//AERR(err);
	
	//     printf("structplusmal %d %d\n", m, cov->nsub);
	// PMI(s->keys[m]);
      }
  
    
      //assert(false);
      return NOERROR;
    }
    
  default :
    SERR2("role '%s' not allowed for '%s'", ROLENAMES[cov->role],
	  NICK(cov));
  }

  return NOERROR;
}


int structplusproc(cov_model *cov, cov_model **newmodel){
  assert(cov->nr == PLUS_PROC);
  return structplusmalproc(cov, newmodel);
}


int structmultproc(cov_model *cov, cov_model **newmodel){
  assert(CovList[cov->nr].check == checkmultproc);
  return structplusmalproc(cov, newmodel);
}

int initplusmalproc(cov_model *cov, gen_storage VARIABLE_IS_NOT_USED *s){
  int i, err,
    vdim = cov->vdim[0];
 assert(cov->vdim[0] == cov->vdim[1]);

  for (i=0; i<vdim; i++)   
    cov->mpp.maxheights[i] = RF_NA;
  if (cov->Splus == NULL) BUG;

  if (cov->role == ROLE_GAUSS) {
 
    for (i=0; i<cov->nsub; i++) {
      //printf("i=%d\n", i);
      cov_model *sub = cov->Splus == NULL ? cov->sub[i] : cov->Splus->keys[i];
      assert(cov->sub[i]->Sgen==NULL);
      cov->sub[i]->Sgen = (gen_storage *) MALLOC(sizeof(gen_storage));
      if ((err = INIT(sub, 0, cov->sub[i]->Sgen)) != NOERROR) {
	return err;
      }
      sub->simu.active = true;
    }
    cov->simu.active = true;
    return NOERROR;
  }
    
  else {
    BUG;
  }

    return ERRORFAILED;

    // assert(false);
}

 
int initplusproc(cov_model *cov, gen_storage VARIABLE_IS_NOT_USED *s){
  int err;
  if ((err = initplusmalproc(cov, s)) != NOERROR) return err;

  if (cov->role == ROLE_GAUSS) {
    cov->fieldreturn = cov->Splus != NULL;
    cov->origrf = false;
    if (cov->Splus != NULL) cov->rf = cov->Splus->keys[0]->rf;
     
    return NOERROR;
  }
     
  else {
    BUG;
  }

  return ERRORFAILED;
}


void doplusproc(cov_model *cov, gen_storage VARIABLE_IS_NOT_USED *s) {
  int m, i,
    total = cov->prevloc->totalpoints * cov->vdim[0];
  double *res = cov->rf;
 assert(cov->vdim[0] == cov->vdim[1]);

  if (cov->role == ROLE_GAUSS && cov->method==SpectralTBM) {
    ERR("error in doplus with spectral");
  }
  assert(cov->Splus != NULL);
  
  for (m=0; m<cov->nsub; m++) {
    cov_model *key = cov->Splus->keys[m],
      *sub = cov->sub[m];
    double *keyrf = key->rf;
    DO(key, sub->Sgen);
    if (m > 0)
      for(i=0; i<total; i++) res[i] += keyrf[i];
  }
  return;
}



#define MULTPROC_COPIES 0
int checkmultproc(cov_model *cov) {
  int err;
  kdefault(cov, MULTPROC_COPIES, GLOBAL.special.multcopies);
  if ((err = checkplusmalproc(cov)) != NOERROR) {
    return err;
  }
  return NOERROR;
}


int initmultproc(cov_model *cov, gen_storage VARIABLE_IS_NOT_USED *s){
  int  err;

  if ((err = initplusmalproc(cov, s)) != NOERROR) {
    BUG;
   return err;
  }

  if (cov->role == ROLE_GAUSS) {
    FieldReturn(cov);
    return NOERROR;
  }
     
  else {
    BUG;
  }

  return ERRORFAILED;
}



void domultproc(cov_model *cov, gen_storage VARIABLE_IS_NOT_USED *s) {
  int m, i,
    total = cov->prevloc->totalpoints * cov->vdim[0];
  double *res = cov->rf;
 assert(cov->vdim[0] == cov->vdim[1]);

  if (cov->role == ROLE_GAUSS && cov->method==SpectralTBM) {
    ERR("error in do_mult with spectral");
  }
  assert(cov->Splus != NULL);

  for(i=0; i<total; res[i++] = 0.0);

  
  for (m=0; m<cov->nsub; m++) {
    cov_model *key = cov->Splus->keys[m],
      *sub = cov->sub[m];
    double *keyrf = key->rf;
    DO(key, sub->Sgen);
      for(i=0; i<total; i++) res[i] += keyrf[i];
  }
  return;
}



void rangemultproc(cov_model VARIABLE_IS_NOT_USED *cov, range_type* range){
  range->min[MULTPROC_COPIES] = 1.0;
  range->max[MULTPROC_COPIES] = RF_INF;
  range->pmin[MULTPROC_COPIES] = 1.0;
  range->pmax[MULTPROC_COPIES] = 1000;
  range->openmin[MULTPROC_COPIES] = false;
  range->openmax[MULTPROC_COPIES] = true;
}


// $power

void PowSstat(double *x, cov_model *cov, double *v){
  logPowSstat(x, cov, v, NULL);

  //printf("powscale %f %f\n", x[0], v);
  
}

void logPowSstat(double *x, cov_model *cov, double *v, double *sign){
  cov_model *next = cov->sub[POW_SUB];
  double 
    factor,
    var = P0(POWVAR),
    scale =P0(POWSCALE), 
    p = P0(POWPOWER),
    invscale = 1.0 / scale;
  int i,
    vdim = cov->vdim[0],
    vdimSq = vdim *vdim,
    xdimown = cov->xdimown;
  assert(cov->vdim[0] == cov->vdim[1]);
  ALLOC_DOLLAR(z, xdimown);

  for (i=0; i < xdimown; i++) z[i] = invscale * x[i];
  if (sign==NULL) {
    COV(z, next, v);
    factor = var * pow(scale, p);
    
    //printf("pows=%f %f %f\n", scale, var, factor);

    for (i=0; i<vdimSq; i++) v[i] *= factor; 
  } else {
    LOGCOV(z, next, v, sign);
    factor = log(var) + p * log(scale);
    for (i=0; i<vdimSq; i++) v[i] += factor; 
  }
}

void PowSnonstat(double *x, double *y, cov_model *cov, double *v){
  logPowSnonstat(x, y, cov, v, NULL);
}

void logPowSnonstat(double *x, double *y, cov_model *cov, double *v, 
		 double *sign){
  cov_model *next = cov->sub[POW_SUB];
  double 
    factor,
    var = P0(POWVAR),
    scale =P0(POWSCALE), 
    p = P0(POWPOWER),
     invscale = 1.0 / scale;
  int i,
    vdim = cov->vdim[0],
    vdimSq = vdim * vdim,
    xdimown = cov->xdimown;
  assert(cov->vdim[0] == cov->vdim[1]);
  ALLOC_DOLLARY(z1, z2, xdimown);

  for (i=0; i<xdimown; i++) {
    z1[i] = invscale * x[i];
    z2[i] = invscale * y[i];
  }

  if (sign==NULL) {
    NONSTATCOV(z1, z2, next, v);
    factor = var * pow(scale, p);
    for (i=0; i<vdimSq; i++) v[i] *= factor; 
  } else {
    LOGNONSTATCOV(z1, z2, next, v, sign);
    factor = log(var) + p * log(scale);
    for (i=0; i<vdimSq; i++) v[i] += factor; 
  }
}

 
void inversePowS(double *x, cov_model *cov, double *v) {
  cov_model *next = cov->sub[POW_SUB];
  int i,
    vdim = cov->vdim[0],
    vdimSq = vdim * vdim;
  double y, 
    scale =P0(POWSCALE),
    p = P0(POWPOWER),
    var = P0(POWVAR);
 assert(cov->vdim[0] == cov->vdim[1]);

  y= *x / (var * pow(scale, p)); // inversion, so variance becomes scale
  if (CovList[next->nr].inverse == ErrCov) BUG;
  INVERSE(&y, next, v);
 
  for (i=0; i<vdimSq; i++) v[i] *= scale; 
}


int TaylorPowS(cov_model *cov) {
  cov_model 
    *next = cov->sub[POW_SUB];
  int i;
  double scale = PisNULL(POWSCALE) ? 1.0 : P0(POWSCALE);
  cov->taylorN = next->taylorN;  
  for (i=0; i<cov->taylorN; i++) {
    cov->taylor[i][TaylorPow] = next->taylor[i][TaylorPow];
    cov->taylor[i][TaylorConst] = next->taylor[i][TaylorConst] *
      P0(POWVAR) * pow(scale, P0(POWPOWER) - next->taylor[i][TaylorPow]);   

    //printf("TC %f %f\n",  cov->taylor[i][TaylorConst], next->taylor[i][TaylorConst]);

  }
  
  cov->tailN = next->tailN;  
  for (i=0; i<cov->tailN; i++) {
    cov->tail[i][TaylorPow] = next->tail[i][TaylorPow];
    cov->tail[i][TaylorExpPow] = next->tail[i][TaylorExpPow];
    cov->tail[i][TaylorConst] = next->tail[i][TaylorConst] *
      P0(POWVAR) * pow(scale, P0(POWPOWER) - next->tail[i][TaylorPow]);   
    cov->tail[i][TaylorExpConst] = next->tail[i][TaylorExpConst] *
      pow(scale, -next->tail[i][TaylorExpPow]);
  }
  return NOERROR;
}

int checkPowS(cov_model *cov) {
  // hier kommt unerwartet  ein scale == nan rein ?!!
  cov_model 
    *next = cov->sub[POW_SUB];
  int err, 
    tsdim = cov->tsdim,
    xdimown = cov->xdimown,
    xdimNeu = xdimown;

  if (!isCartesian(cov->isoown)) return ERRORNOTCARTESIAN;
    
  kdefault(cov, POWVAR, 1.0);
  kdefault(cov, POWSCALE, 1.0);
  kdefault(cov, POWPOWER, 0.0);
  if ((err = checkkappas(cov)) != NOERROR) {
    return err;
  }
  
  if ((err = CHECK(next, tsdim, xdimNeu, cov->typus, cov->domown,
		   cov->isoown, SUBMODEL_DEP, cov->role)) != NOERROR) {
    return err;
  }

  setbackward(cov, next);
  if ((err = TaylorPowS(cov)) != NOERROR) return err;

  DOLLAR_STORAGE;
  
  return NOERROR;
}



bool TypePowS(Types required, cov_model *cov, int depth) {
  cov_model *next = cov->sub[0];

  // printf("\n\n\nxxxxxxxxxxxx %d %s %s\n\n", TypeConsistency(required, next), TYPENAMES[required], NICK(next));
  
  // if (required == ProcessType) crash(); //assert(false);

  return (required==TcfType || required==PosDefType || required==NegDefType
	  || required==ShapeType || required==TrendType 
	  || required==ProcessType || required==GaussMethodType)
    && TypeConsistency(required, next, depth-1);
}


void rangePowS(cov_model *cov, range_type* range){
  range->min[POWVAR] = 0.0;
  range->max[POWVAR] = RF_INF;
  range->pmin[POWVAR] = 0.0;
  range->pmax[POWVAR] = 100000;
  range->openmin[POWVAR] = false;
  range->openmax[POWVAR] = true;

  range->min[POWSCALE] = 0.0;
  range->max[POWSCALE] = RF_INF;
  range->pmin[POWSCALE] = 0.0001;
  range->pmax[POWSCALE] = 10000;
  range->openmin[POWSCALE] = true;
  range->openmax[POWSCALE] = true;

  range->min[POWPOWER] = RF_NEGINF;
  range->max[POWPOWER] = RF_INF;
  range->pmin[POWPOWER] = -cov->tsdim;
  range->pmax[POWPOWER] = +cov->tsdim;
  range->openmin[POWPOWER] = true;
  range->openmax[POWPOWER] = true;
 }



void PowScaleToLoc(cov_model *to, cov_model *from, int VARIABLE_IS_NOT_USED depth) {
  assert(!PARAMisNULL(to, LOC_SCALE) && !PARAMisNULL(from, POWSCALE));
  PARAM(to, LOC_SCALE)[0] = PARAM0(from, POWSCALE);
}

int structPowS(cov_model *cov, cov_model **newmodel) {
  cov_model
    *next = cov->sub[POW_SUB],
    *scale = cov->kappasub[POWSCALE];
  int err; 
  // cov_model *next;



  //   printf("%s %d %s\n", ROLENAMES[cov->role], cov->role, TYPENAMES[cov->typus]);
  //  if (cov->role != ROLE_GAUSS) {
    //  APMI(cov->calling);
    //crash();
  //  }
  if (!next->deterministic) SERR("random shapes not programmed yet");

  switch (cov->role) {
  case ROLE_SMITH :  case ROLE_GAUSS :
    ASSERT_NEWMODEL_NOT_NULL;
    
    if ((err = STRUCT(next, newmodel)) > NOERROR) return err;
    
    addModel(newmodel, POWER_DOLLAR);
    if (!PisNULL(POWVAR)) kdefault(*newmodel, POWVAR, P0(POWVAR));
    if (!PisNULL(POWSCALE)) kdefault(*newmodel, POWSCALE, P0(POWSCALE));
    if (!PisNULL(POWPOWER)) kdefault(*newmodel, POWPOWER, P0(POWPOWER));
    
    break;
  case ROLE_MAXSTABLE : {
    ASSERT_NEWMODEL_NOT_NULL;
    
    if ((err = STRUCT(next, newmodel)) > NOERROR) return err;
    
    if (!isRandom(scale)) SERR("unstationary scale not allowed yet");
    addModel(newmodel, LOC);
    addSetDistr(newmodel, scale, PowScaleToLoc, true, MAXINT);
  }
    break;
  default :
    //  PMI(cov, "structS");
    SERR2("'%s': changes in scale/variance not programmed yet for '%s'", 
	  NICK(cov), ROLENAMES[cov->role]);      
  }  
   
  return NOERROR;
}




int initPowS(cov_model *cov, gen_storage *s){
  // am liebsten wuerde ich hier die Koordinaten transformieren;
  // zu grosser Nachteil ist dass GetDiameter nach trafo 
  // grid def nicht mehr ausnutzen kann -- umgehbar?!
  cov_model *next = cov->sub[POW_SUB],
    *varM = cov->kappasub[POWVAR],
    *scaleM = cov->kappasub[POWSCALE];
  int 
    vdim = cov->vdim[0],
    nm = cov->mpp.moments,
    nmvdim = (nm + 1) * vdim,
    err = NOERROR;
  bool 
    maxstable = hasExactMaxStableRole(cov);// Realisationsweise 
  assert(cov->vdim[0] == cov->vdim[1]);


  assert(cov->key == NULL || ({PMI(cov);false;}));//   
  
  if (hasAnyShapeRole(cov)) { // !! ohne maxstable selbst !!
    double
      var[MAXMPPVDIM],  
      p = P0(POWPOWER),
      scale = P0(POWSCALE);
    int i,
      intp = (int) p,
      dim = cov->tsdim;

    

    // Achtung I-NIT_RANDOM ueberschreibt mpp.* !!
    if (varM != NULL) {
      int nm_neu = nm == 0 && !maxstable ? 1 : nm;
      if ((err = INIT_RANDOM(varM, nm_neu, s, P(POWVAR))) != NOERROR) 
	return err;
      int nmP1 = varM->mpp.moments + 1;
      for (i=0; i<vdim; i++) {
	int idx = i * nmP1;
	var[i] = maxstable ? P0(DVAR) : varM->mpp.mM[idx + 1];      
      }
    } else for (i=0; i<vdim; var[i++] = P0(POWVAR));

    if (scaleM != NULL) {
      if (p != intp)
	SERR1("random scale can be initialised only for integer values of '%s'",
	     KNAME(POWPOWER));
      if (dim + intp < 0) SERR("negative power cannot be calculated yet");
      int idx_s = maxstable ? nm : dim + nm + intp < 1 ? 1 : dim + nm + intp;
      if ((err = INIT_RANDOM(scaleM, idx_s, s, P(POWSCALE))) != NOERROR) return err;
      assert(scaleM->mpp.moments == 1);
      scale = maxstable ? P0(DSCALE) : scaleM->mpp.mM[1];      
    }
    if ((err = INIT(next, nm, s)) != NOERROR) return err;


    for (i=0; i < nmvdim; i++) {
      cov->mpp.mM[i] = next->mpp.mM[i];
      cov->mpp.mMplus[i] = next->mpp.mMplus[i];
      //printf("%s I=%d %f %f\n", NICK(cov->calling), i,  cov->mpp.mMplus[i], cov->mpp.mM[i]);
   }


    if (varM != NULL && !maxstable) {
      int j,
	nmP1 = varM->mpp.moments + 1;
      for (j=0; j<vdim; j++) {
	int idx = j * nmP1;
	for (i=0; i <= nm; i++) {
	  cov->mpp.mM[i] *= varM->mpp.mM[idx + i];
	  cov->mpp.mMplus[i] *= varM->mpp.mMplus[idx + i];
	}
      }
    } else {
      int j, k;
      double pow_var;
      for (k=j=0; j<vdim; j++) { 
	pow_var = 1.0;
	for (i=0; i<=nm; i++, pow_var *= var[j], k++) {	
	  cov->mpp.mM[k] *= pow_var;
	  cov->mpp.mMplus[k] *= pow_var;
	  //printf("%s i=%d %f p=%f %f\n", NICK(cov->calling), i, var, pow_var, cov->mpp.mM[i]);
	}
      }
    }

    if (scaleM != NULL && !maxstable) {
      if (dim + nm * intp < 0 || dim + intp * nm > scaleM->mpp.moments) 
	SERR("moments cannot be calculated");
      assert(scaleM->vdim[0] == 1 && scaleM->vdim[1] == 1 );
      for (i=0; i <= nm; i++) {
	int idx = dim + i * intp;
	cov->mpp.mM[i] *= scaleM->mpp.mM[idx];
	cov->mpp.mMplus[i] *= scaleM->mpp.mMplus[idx];
      }
    } else {
      int j,k;
      double 
	pow_scale,
	pow_s = pow(scale, dim),
	pow_p = pow(scale, p);
      for (k=j=0; j<vdim; j++) { 
	pow_scale = pow_s;
	for (i=0; i <= nm; i++, pow_scale *= pow_p, k++) {
	  cov->mpp.mM[k] *= pow_scale;
	  cov->mpp.mMplus[k] *= pow_scale;
	}
    	//printf("%s i=%d s=%e %e\n", NICK(cov->calling), i, pow_scale, cov->mpp.mM[i]);
      }
    }
 
    if (R_FINITE(cov->mpp.unnormedmass)) {
      if (vdim > 1) BUG;
      cov->mpp.unnormedmass = next->mpp.unnormedmass * var[0] / pow(scale, p);
    } else {
      double pp = pow(scale, -p);
      for (i=0; i<vdim; i++)   
	cov->mpp.maxheights[i] = next->mpp.maxheights[i] * var[i] * pp;
    }
 
    //    printf("ini pows %f %f %d %f\n",   
    //	   cov->mpp.maxheight , next->mpp.maxheight ,
    //	   varM == NULL, varM == NULL ? P0(POWVAR) : 1.0);
 
    //    cov->mpp.refradius *= scale;
    //cov->mpp.refsd *= scale;
  }


  else if (cov->role == ROLE_GAUSS) {  
    if ((err=INIT(next, 0, s)) != NOERROR) return err;
  }

  else if (cov->role == ROLE_BASE) {
    if ((err=INIT(next, 0, s)) != NOERROR) return err;
    
  }

  else SERR("initiation of scale and/or variance failed");

 
  if ((err = TaylorPowS(cov)) != NOERROR) return err;

  return NOERROR;
}


void doPowS(cov_model *cov, gen_storage *s){
 
  if (hasAnyShapeRole(cov)) {
    cov_model *next = cov->sub[POW_SUB];
         
    DO(next, s);// nicht gatternr
    double factor = P0(POWVAR) / pow(P0(POWSCALE), P0(POWPOWER));
    int i,
      vdim = cov->vdim[0];
    assert(cov->vdim[0] == cov->vdim[1]);
    for (i=0; i<vdim; i++)   
      cov->mpp.maxheights[i] = next->mpp.maxheights[i] * factor;
    return;
  }
  
  BUG;
}