RFtbm.cc
/*
Authors
Martin Schlather, Martin.Schlather@uni-bayreuth.de
Simulation of a random field by turning bands;
see RFspectral.cc for spectral turning bands
Copyright (C) 2001 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 2
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 <stdlib.h>
#include <sys/timeb.h>
#include <unistd.h>
#include <assert.h>
#include "RFsimu.h"
#include "RFCovFcts.h"
#define MAXNN 100000000.0 /* max number of points simulated on a line */
int TBM2_LINES=60;// number of lines simulated
int TBM3D2_LINES=500;
int TBM3D3_LINES=500;
Real TBM2_LINESIMUFACTOR=2.0; /*factor by which simulation on the line is finer
than the grid */
Real TBM2_LINESIMUSTEP=0.0; /* grid lag defined absolutely */
Real TBM3D2_LINESIMUFACTOR=2.0;
Real TBM3D2_LINESIMUSTEP=0.0;
Real TBM3D3_LINESIMUSTEP=0.0;
Real TBM3D3_LINESIMUFACTOR=2.0;
bool TBMCE_FORCE=false; /* if TRIALS times failed then do as if it was OK?
see also CIRCEMBED_* */
int TBMCE_TRIALS=3;/* how often size is doubled in case of failure before
giving up */
int TBMCE_MTOT=0; /* minimal mtot */
#ifdef doubleReal
Real TBMCE_TOL_RE=-1e-7;
Real TBMCE_TOL_IM=1e-3;
#else
Real TBMCE_TOL_RE=-1e-5;
Real TBMCE_TOL_IM=1e-3;
#endif
SimulationType TBM_METHOD=Nothing;
typedef struct TBM_storage {
simu1dim method; /* simulation method on the line (up to now only
circulantembedding) */
int nn;
Real linesimuscale,half[MAXDIM];
int mtot;
double *c,*d,*simuline;
FFT_storage FFT;
} TBM_storage;
void TBM_destruct(void **S)
{
if (*S!=NULL) {
TBM_storage *x;
x = *((TBM_storage**)S);
if (x->c!=NULL) free(x->c);
if (x->d!=NULL) free(x->d);
if (x->simuline!=NULL) free(x->simuline);
FFT_destruct(&(x->FFT));
free(*S);
*S = NULL;
}
}
// the following two procedures are the same as in RFcircembed, except simplified
// to dimension 1
void do_circulantembedding1dim(key_type *key,Real *res END_WITH_GSLRNG)
{
int dim, i, kk,
index,halfmtot,zeroORmiddle,
modtot,mtot;
Real *d,*c,UU,VV,XX,YY,invsqrtmtot;
int TwoIndex,TwoKK,twoRealmtot,nn;
// int OneDimMaxm=67108864;
bool first;
#ifdef RF_GSL
assert(RANDOM!=NULL);
#endif
assert(key->active);
nn = ((TBM_storage*)key->S)->nn;
mtot = ((TBM_storage*)key->S)->mtot;
// PRINTF(" %d mtot =%d \n",nn,mtot);
c= ((TBM_storage*)key->S)->c;
d= ((TBM_storage*)key->S)->d;
twoRealmtot = 2*sizeof(Real) * mtot;
halfmtot = mtot/2;
modtot = halfmtot + 1;
invsqrtmtot = 1/sqrt((Real) mtot);
dim=1;
/* now the Gaussian r.v. have to defined and multiplied with sqrt(FFT(c))*/
for (index=0; index<modtot; index++) {
/* about half of the loop is clearly unnecessary,
as pairs of r.v. are created, see inequality * below;
therefore modtot is used instead of mtot;
as the procedure is complicated, it is easier to let the loop run
over a bit more values of i than using exactly mtot/2 values */
zeroORmiddle = 1;
kk=0; /* kk takes the role of k[l] in Prop. 3; but here kk is already
transformed from matrix indices to a vector index */
if (index==0 || index==halfmtot) {kk += index;}
else {
zeroORmiddle=0;
kk += mtot-index;
}
if (zeroORmiddle){ /* step 2 of Prop. 3 */
if (index<halfmtot) {
/* create two Gaussian random variables XX & YY with variance 1/2 */
UU= sqrt(-log(1.0 - UNIFORM_RANDOM));
VV = TWOPI * UNIFORM_RANDOM;
XX = UU * sin(VV);
YY = UU * cos(VV);
kk = index + halfmtot;
TwoIndex = index <<1;
d[TwoIndex] = c[TwoIndex] * XX; d[TwoIndex+1]=0.0;
TwoKK=kk <<1;
d[TwoKK] = c[TwoKK] * YY; d[TwoKK+1]=0.0;
}
} else { /* step 3 of Prop. 3 */
UU= sqrt(-log(1.0 - UNIFORM_RANDOM));
VV = TWOPI * UNIFORM_RANDOM;
XX = UU * sin(VV);
YY = UU * cos(VV);
TwoIndex= index <<1;
d[TwoIndex+1]=(d[TwoIndex]=c[TwoIndex])*YY; d[TwoIndex]*=XX;
TwoKK=kk <<1;
d[TwoKK+1]=(d[TwoKK]=c[TwoKK]) * (-YY); d[TwoKK]*=XX;
}
}
fastfourier(d,&mtot,dim,FALSE,&(((TBM_storage*)key->S)->FFT));
first =true;
for (kk=i=0;i<nn;i++,kk++){
res[i]= (Real) d[kk++] *invsqrtmtot;
if ((fabs(d[kk])>TBMCE_TOL_IM) && (first)){
if (GENERAL_PRINTLEVEL>=2)
PRINTF("IMAGINARY PART <> 0; %d %f \n",kk,d[kk]); first=false;
}
}
}
int init_circulantembedding1dim(key_type *key, Real *p)
{
int dim=1, i,j,index,halfmtot, mtot,TwoI,trials,nn;
Real *c;
int twoRealmtot,error;
int MAXM=67108864;
bool positivedefinite;
covfct cov;
c = NULL;
nn = ((TBM_storage*)key->S)->nn;
mtot = 1 << (1 + (int) ceil(log(((Real) nn)-EPSILON1000) * INVLOG2)); // nn<=2^(mtot-1)
if (TBMCE_MTOT>mtot) {mtot=TBMCE_MTOT;}
if (mtot>MAXM) { error=ERRORFAILED; goto ErrorHandling;}
if (key->method==TBM2) { cov=key->cov->cov_tbm2;}
else {
assert(key->method==TBM3);
cov=key->cov->cov_tbm3;
}
positivedefinite = false;
trials =0;
while (!positivedefinite && (trials<TBMCE_TRIALS)){
halfmtot = mtot/2;
index=1-halfmtot;
trials++;
if (GENERAL_PRINTLEVEL>=2) {
PRINTF("trial %d of %d\n",trials,TBMCE_TRIALS);
PRINTF("mtot=%d\n",mtot);
}
twoRealmtot = 2*sizeof(Real) * mtot;
if ((c =(Real*) malloc(twoRealmtot)) ==0){
error=ERRORMEMORYALLOCATION; goto ErrorHandling;
}
for (i=0;i<mtot;i++){ /* fill in c(h) column-wise*/
j= 2 * ((index+mtot) % mtot);
c[j]=cov((Real) index,p); c[j+1]=0;
if (++index>halfmtot) {index=1-halfmtot;}
}
if (error=(fastfourier(c,&mtot,(int) dim,TRUE,&(((TBM_storage*)key->S)->FFT))))
goto ErrorHandling;
/* check if positive definite. If not enlarge and restart */
if (!TBMCE_FORCE || (trials<TBMCE_TRIALS)) {
i=0; TwoI = 0;
while ((i<mtot)&&(positivedefinite=(c[TwoI]>=TBMCE_TOL_RE &&
fabs(c[TwoI+1])<TBMCE_TOL_IM)))
{// Approximation if in tolerance bounds
i++; TwoI += 2;
}
if (!positivedefinite) {
if (GENERAL_PRINTLEVEL>=2)
PRINTF(" nonpos %d %f %f \n",i,c[TwoI],c[TwoI+1]);
free(c); c=NULL;
mtot <<= 1;
if (mtot>MAXM) { error=ERRORFAILED; goto ErrorHandling;}
}
} else {if (GENERAL_PRINTLEVEL>=2) PRINTF("forced\n");}
}
// note: mtot<<=1 above; so now mtot equals number of allocated units of d
if (positivedefinite || TBMCE_FORCE) {
for(i=0,TwoI=0;i<mtot;i++,TwoI+=2) {
c[TwoI] = (c[TwoI]>0.0 ? sqrt(c[TwoI]) : 0.0);
c[TwoI+1] = 0.0;
}
} else {
error=ERRORFAILED; goto ErrorHandling;
}
if ((((TBM_storage*)key->S)->simuline= (Real *)malloc(nn * sizeof(Real)))==0){
error=ERRORMEMORYALLOCATION; goto ErrorHandling;
}
if ((((TBM_storage*)key->S)->d = (Real *)malloc(twoRealmtot)) ==0){
error=ERRORMEMORYALLOCATION; goto ErrorHandling;
}
((TBM_storage*)key->S)->mtot = mtot;
((TBM_storage*)key->S)->c = c;
return 0;
ErrorHandling:
if (c!=NULL) free(c);
assert(key->destruct!=NULL);
key->destruct(&key->S);
key->destruct=NULL;
return error;
}
void SetParamTBMCE(int *action,int *force,Real *tolRe,Real *tolIm,int *trials,
int *mtot) {
switch (*action) {
case 0 :
TBMCE_FORCE=(bool) *force;
TBMCE_TRIALS = *trials;
if (TBMCE_TRIALS<1) {
TBMCE_TRIALS=1;
if (GENERAL_PRINTLEVEL>0)
PRINTF("\nWARNING! TBMCE_TRIALS had been less than 1\n");
}
TBMCE_TOL_RE=*tolRe;
if (TBMCE_TOL_RE>0) {
TBMCE_TOL_RE=0;
if (GENERAL_PRINTLEVEL>0)
PRINTF("\nWARNING! TBMCE_TOL_RE had been positive.\n");
}
TBMCE_TOL_IM=*tolIm;
if (TBMCE_TOL_IM<0) {
TBMCE_TOL_IM=0;
if (GENERAL_PRINTLEVEL>0)
PRINTF("\nWARNING! TBMCE_TOL_IM had been neagtive.\n");
}
TBMCE_MTOT=*mtot;
if (TBMCE_MTOT>0) {
TBMCE_MTOT = 1 << (1+(int) (log(((double) TBMCE_MTOT)-0.1)/log(2.0)));
if ((TBMCE_MTOT!=*mtot) && (GENERAL_PRINTLEVEL>0))
PRINTF("\nWARNING! TBMCE_MTOT set to %d.\n",TBMCE_MTOT);
}
break;
case 1 :
*force=TBMCE_FORCE;
*tolRe=TBMCE_TOL_RE;
*tolIm=TBMCE_TOL_IM;
*trials=TBMCE_TRIALS;
*mtot=TBMCE_MTOT;
if (GetNotPrint) break;
case 2 :
PRINTF("\nTBM_CE\n======\nforce=%d\ntol_Re=%e\ntol_Im=%e\ntrials=%d\nmtot=%d\n",
TBMCE_FORCE,TBMCE_TOL_RE,TBMCE_TOL_IM,TBMCE_TRIALS,TBMCE_MTOT);
break;
default : PRINTF(" unknown action\n");
}
}
void SetParamTBM2(int *action,int *nLines, Real *linesimufactor,
Real *linesimustep) {
switch (*action) {
case 0 :
TBM2_LINES=*nLines;
if ((*linesimufactor>0.0) ^ (*linesimustep>0.0)) {
TBM2_LINESIMUFACTOR=*linesimufactor;
TBM2_LINESIMUSTEP=*linesimustep;
} else {
PRINTF(" Exactly one of `linesimufactor' and `linesimustep' should be positive.");
}
break;
case 1 :
*nLines=TBM2_LINES;
*linesimufactor=TBM2_LINESIMUFACTOR;
*linesimustep=TBM2_LINESIMUSTEP;
if (GetNotPrint) break;
case 2:
PRINTF("\nTBM2\n====\nLines=%d\nlinesimufactor=%f\nlinesimustep=%f\n",
TBM2_LINES,TBM2_LINESIMUFACTOR,TBM2_LINESIMUSTEP);
break;
default : PRINTF(" unknown action\n");
}
}
void SetParamTBM3D2(int *action,int *nLines, Real *linesimufactor,
Real *linesimustep){
switch (*action) {
case 0 :
TBM3D2_LINES=*nLines;
if ((*linesimufactor>0.0)^(*linesimustep>0.0)) {
TBM3D2_LINESIMUFACTOR=*linesimufactor;
TBM3D2_LINESIMUSTEP=*linesimustep;
} else {
PRINTF(" Exactly one of `linesimufactor' and `linesimustep' should be positive");
}
break;
case 1 :
*nLines=TBM3D2_LINES;
*linesimufactor=TBM3D2_LINESIMUFACTOR;
*linesimustep=TBM3D2_LINESIMUSTEP;
if (GetNotPrint) break;
case 2:
PRINTF("\nTBM3D2\n======\nLines=%d\nlinesimufactor=%f\nlinesimustep=%f\n",
TBM3D2_LINES,TBM3D2_LINESIMUFACTOR,TBM3D2_LINESIMUSTEP);
break;
default : PRINTF(" unknown action\n");
}
}
void SetParamTBM3D3(int *action,int *nLines, Real *linesimufactor,
Real *linesimustep){
switch (*action) {
case 0 :
TBM3D3_LINES=*nLines;
if ((*linesimufactor>0.0)^(*linesimustep>0.0)) {
TBM3D3_LINESIMUFACTOR=*linesimufactor;
TBM3D3_LINESIMUSTEP=*linesimustep;
} else {
PRINTF(" Exactly one of `linesimufactor' and `linesimustep' should be positive");
}
break;
case 1 :
*nLines=TBM3D3_LINES;
*linesimufactor=TBM3D3_LINESIMUFACTOR;
*linesimustep=TBM3D3_LINESIMUSTEP;
if (GetNotPrint) break;
case 2:
PRINTF("\nTBM3D3\n======\nLines=%d\nlinesimufactor=%f\nlinesimustep=%f\n",
TBM3D3_LINES,TBM3D3_LINESIMUFACTOR,TBM3D3_LINESIMUSTEP);
break;
default : PRINTF(" unknown action\n");
}
}
void SetParamTBM(int *action, int *tbm_method) {
switch (*action) {
case 0 :
if ((*tbm_method<=(int) Nothing) && (*tbm_method>=0)) {
SimulationType m;
m=(SimulationType) *tbm_method;
if ((m!=TBM2)&&(m!=TBM3)&&(m!=SpectralTBM)){
TBM_METHOD = m;
return;
}
}
if (GENERAL_PRINTLEVEL>0) {
PRINTF("Specified method [%d] is not allowed. Method set to `Nothing'\n",
*tbm_method);
}
TBM_METHOD = Nothing;
break;
case 1 :
*tbm_method = (int) TBM_METHOD;
if (GetNotPrint) break;
case 2 :
PRINTF("\nTBM:\n====\nuser defined method=%d\n",TBM_METHOD);
break;
default : PRINTF(" unknown action\n");
}
}
int init_turningbands(key_type *key)
{
int error,i,d;
Real p[TOTAL_PARAM],linesimufactor,linesimustep,dist,mindelta;
assert((key->covnr>=0) &&(key->covnr<currentNrCov));
assert(key->active==false);
assert(key->param[VARIANCE]>=0.0);
assert(key->param[SILL]>=key->param[VARIANCE]);
assert( fabs(key->param[SCALE]*key->param[INVSCALE]-1.0) < EPSILON);
assert((key->S==NULL) && (key->destruct==NULL));
if(key->cov->check!=NULL) if(error=key->cov->check(key)) goto ErrorHandling;
if ((key->S=malloc(sizeof(TBM_storage)))==0){
error=ERRORMEMORYALLOCATION; goto ErrorHandling;
}
((TBM_storage*)key->S)->c =NULL;
((TBM_storage*)key->S)->d =NULL;
((TBM_storage*)key->S)->simuline=NULL;
FFT_NULL(&(((TBM_storage*)key->S)->FFT));
key->destruct = TBM_destruct;
if (key->dim==2) {
if (key->method==TBM2) {
linesimufactor=TBM2_LINESIMUFACTOR;
linesimustep = TBM2_LINESIMUSTEP;
}
else {
assert(key->method==TBM3);
linesimufactor=TBM3D2_LINESIMUFACTOR;
linesimustep = TBM3D2_LINESIMUSTEP;
}
} else if (key->dim==3) {
assert(key->method==TBM3);
linesimufactor=TBM3D3_LINESIMUFACTOR;
linesimustep = TBM3D3_LINESIMUSTEP;
} else assert(false);
for (i=0;i<TOTAL_PARAM;i++) {p[i]=key->param[i];}
if (key->grid) {
{
register Real dummy;
dist = 0.0;
for (d=0; d<key->dim; d++) {
dummy = key->x[d][XSTEP] * (Real) (key->length[d] - 1);
dist += dummy * dummy;
}
}
if (linesimustep>0.0) {
((TBM_storage*)key->S)->linesimuscale = 1/linesimustep;
}
else {
mindelta = RF_INF;
for (d=0; d<key->dim; d++) {
if ((key->x[d][XSTEP]<mindelta) && (key->x[d][XSTEP]>0))
{mindelta=key->x[d][XSTEP];}
}
((TBM_storage*)key->S)->linesimuscale = linesimufactor/mindelta;
}
} else { // not a grid
Real min[MAXDIM],max[MAXDIM];
int dim,j,d;
dim = key->dim;
if (linesimustep>0.0) {
((TBM_storage*)key->S)->linesimuscale = 1/linesimustep;
}
else {
mindelta = RF_INF;
for (i=0;i<key->totallength;i++) {
for (j=0;j<i; j++) {
register Real diff,dist;
for(dist=0.0,d=0;d<dim;d++) {
diff = key->x[d][i]- key->x[d][j];
dist += diff * diff;
}
if (dist<mindelta) mindelta=dist;
}
}
((TBM_storage*)key->S)->linesimuscale = linesimufactor/sqrt(mindelta);
}
for (d=0;d<MAXDIM;d++) {min[d]=RF_INF; max[d]=RF_NEGINF;}
for (i=0;i<key->totallength;i++) {
for (d=0;d<dim;d++) {
if (key->x[d][i]<min[d]) min[d] = key->x[d][i];
if (key->x[d][i]>max[d]) max[d] = key->x[d][i];
}
}
dist = 0.0;
for (d=0; d<key->dim; d++) {
dist += (max[d]-min[d])*(max[d]-min[d]);
((TBM_storage*)key->S)->half[d] = (max[d]+min[d])/2.0;
}
}
{
Real dummy;
dummy = (sqrt(dist)* ((TBM_storage*)key->S)->linesimuscale);
if (dummy > MAXNN) {error=ERRORNN; goto ErrorHandling;}
((TBM_storage*)key->S)->nn= 5+ (long) dummy;//+5 for safety
}
p[SCALE]=key->param[SCALE] * ((TBM_storage*)key->S)->linesimuscale ;
p[INVSCALE] = 1.0/(p[SCALE]);
if (GENERAL_PRINTLEVEL>=4) {
PRINTF(" nn =%d\n", ((TBM_storage*)key->S)->nn);
PRINTF("TBM invscale= %f; %f %f %f %f \n\n",
p[INVSCALE],linesimustep,key->param[SCALE],
linesimufactor,key->x[0][XSTEP]);
}
// nn=number of pixels to be simulated on the lines [==(length of diagonal
// through grid that is to be simulated) * (linesimufactor=="precision")]
// 5 is for safety
// maybe also loop and run through preference lists ...
switch (key->tbm_method) {
case CircEmbed :
((TBM_storage*)key->S)->method = do_circulantembedding1dim;
error=init_circulantembedding1dim(key,p);
break;
default : error=ERRORNOTPROGRAMMED;
}
if (!error) return 0;
ErrorHandling:
if (key->destruct!=NULL) key->destruct(&key->S);
key->destruct=NULL;
key->active = false;
return error;
}
void do_turningbands(key_type *key,Real *res END_WITH_GSLRNG)
{
Real phi,InvSqrtNlines;
Real *simuline;
long n,nn,nnhalf;
simu1dim simulate;
#ifdef RF_GSL
assert(RANDOM!=NULL);
#endif
assert(key->active);
simulate = ((TBM_storage*)key->S)->method;
nn = ((TBM_storage*)key->S)->nn;
// printf("nnn = %d\n",nn);
nnhalf = nn/2;
simuline = ((TBM_storage*)key->S)->simuline;
{
register long i,endfor;
endfor = key->totallength;
for (i=0;i<endfor;i++) { res[i]=0.0; }
}
if (key->grid) {
Real deltax, yoffset,deltay, nxhalf, nyhalf,xoffset,
linesimuscaleX,linesimuscaleY,linesimuscaleZ;
int nx,zaehler,ny;
nxhalf=((Real)(key->length[0]-1))/2.0;
nyhalf=((Real)(key->length[1]-1))/2.0;
linesimuscaleX = ((TBM_storage*)key->S)->linesimuscale * key->x[0][XSTEP];
linesimuscaleY = ((TBM_storage*)key->S)->linesimuscale * key->x[1][XSTEP];
if (key->method==TBM2) {
Real deltaphi;
InvSqrtNlines=1.0 / sqrt((Real) TBM2_LINES);
deltaphi = PI / (Real) TBM2_LINES;
phi = deltaphi * UNIFORM_RANDOM;
for (n=0;n<TBM2_LINES;n++){
deltax=sin(phi) * linesimuscaleX;
deltay=cos(phi) * linesimuscaleY;
simulate(key,simuline END_WITH_RANDOM); /* depending on deltax and deltay,
shorter lines might be
simulated! -> improvement in
speed?! */
/* then nnhalf must be kept variable, too: */
yoffset= (double)nnhalf - nyhalf * deltay - nxhalf * deltax;
zaehler = 0;
for (ny=0;ny<key->length[1]; ny++) {
xoffset = yoffset;
for(nx=0;nx<key->length[0]; nx++) {
assert((xoffset<nn) && (xoffset>=0) );
res[zaehler++] += simuline[(long) xoffset];
xoffset += deltax;
}
yoffset += deltay;
}
phi += deltaphi;
}
} else {
assert(key->method==TBM3);
if (key->dim==2) { // TBM3
InvSqrtNlines=1.0 / sqrt((Real) TBM3D2_LINES);
for (n=0;n<TBM3D2_LINES;n++){
if ((GENERAL_PRINTLEVEL>=3) && (n % 100 == 0)) PRINTF("%d \n",n);
deltax= UNIFORM_RANDOM;// see Martin's tech rep for details
deltay=
sqrt(1.0-deltax*deltax) * sin(UNIFORM_RANDOM*TWOPI) * linesimuscaleY;
deltax *= linesimuscaleX;
simulate(key,simuline END_WITH_RANDOM);
/* depending on deltax and deltay, shorter lines might be
simulated! -> improvement in speed!
*/
yoffset= (double)nnhalf - nyhalf * deltay - nxhalf * deltax ;
zaehler = 0;
for (ny=0;ny<key->length[1];ny++) {
xoffset = yoffset;
for(nx=0;nx<key->length[0];nx++) {
assert((xoffset<nn) && (xoffset>=0) );
res[zaehler++] += simuline[(long) xoffset];
if (!(fabs(simuline[(long) xoffset])<10000.0))
PRINTF("s:%f ",simuline[(long) xoffset]);
xoffset += deltax;
}
yoffset += deltay;
}
}
} else {
Real dummy,nzhalf,zoffset,deltaz;
int nz;
assert(key->dim==3);
linesimuscaleZ= ((TBM_storage*)key->S)->linesimuscale * key->x[2][XSTEP];
nzhalf=((Real)(key->length[2]-1))/2.0;
InvSqrtNlines=1.0 / sqrt((Real) TBM3D3_LINES);
for (n=0;n<TBM3D3_LINES;n++){
if ((GENERAL_PRINTLEVEL>=3) && (n % 100 == 0)) PRINTF("%d \n",n);
dummy = UNIFORM_RANDOM * PI;
deltaz = sin(dummy) * linesimuscaleZ;
dummy = cos(dummy);
deltay= UNIFORM_RANDOM * TWOPI;
deltax= cos(deltay) * dummy * linesimuscaleX;
deltay= sin(deltay) * dummy * linesimuscaleY;
simulate(key,simuline END_WITH_RANDOM);
zoffset=
(double)nnhalf - nzhalf * deltaz - nyhalf * deltay - nxhalf * deltax;
zaehler = 0;
for (nz=0;nz<key->length[2];nz++) {
yoffset = zoffset;
for (ny=0;ny<key->length[1];ny++) {
xoffset = yoffset;
for(nx=0;nx<key->length[0];nx++) {
///
if ((xoffset>=nn) || (xoffset<0))
PRINTF(" %d %d %d | %d %d %d \n %f %f %f\n %d %f %f %f\n %d %f\n",
nx,ny,nz,key->length[0],key->length[1],key->length[2],
deltax,deltay,deltaz,
nnhalf,nxhalf,nyhalf,nzhalf,
nn,xoffset);
///
assert( (xoffset<nn) && (xoffset>=0) );
res[zaehler++] += simuline[(long) xoffset];
if (!(fabs(simuline[(long) xoffset])<10000.0))
PRINTF("s:%f ",simuline[(long) xoffset]);
xoffset += deltax;
}
yoffset += deltay;
}
zoffset += deltaz;
}
}
}
}
} else { // not grid
Real halfx,halfy,ex,ey,offset,linescale;
halfx = ((TBM_storage*)key->S)->half[0];
halfy = ((TBM_storage*)key->S)->half[1];
linescale = ((TBM_storage*)key->S)->linesimuscale;
int i;
if (key->method==TBM2) {
Real deltaphi;
InvSqrtNlines=1.0 / sqrt((Real) TBM2_LINES);
deltaphi = PI / (Real) TBM2_LINES;
phi = deltaphi * UNIFORM_RANDOM;
for (n=0;n<TBM2_LINES;n++){
ex=sin(phi) * linescale;
ey=cos(phi) * linescale;
simulate(key,simuline END_WITH_RANDOM);
/* depending on deltax and deltay, shorter lines might be simulated!
-> improvement in speed! */
offset= nnhalf - halfy * ey - halfx * ex;
for (i=0;i<key->totallength;i++) {
register long index;
index = (long) (offset + key->x[0][i] * ex + key->x[1][i] * ey);
assert((index<nn) && (index>=0));
res[i] += simuline[index];
}
phi += deltaphi;
}
} else { // TBM3
assert(key->method==TBM3);
if (key->dim==2) {
Real linesimusq;
linesimusq = linescale * linescale;
InvSqrtNlines=1.0 / sqrt((Real) TBM3D2_LINES);
for (n=0;n<TBM3D2_LINES;n++){
if ((GENERAL_PRINTLEVEL>=3) && (n % 100 == 0)) PRINTF("%d \n",n);
ex= UNIFORM_RANDOM * linescale;// see Martin's tech rep for details
ey= sqrt(linesimusq - ex * ex) * sin(UNIFORM_RANDOM * TWOPI);
simulate(key,simuline END_WITH_RANDOM);
/* depending on ex and ey, shorter lines might be simulated! ->
improvement in speed! */
offset= (double) nnhalf - halfy * ey - halfx * ex;
for (i=0;i<key->totallength;i++) {
register long index;
index = (long) (offset + key->x[0][i] * ex + key->x[1][i] * ey);
assert((index<nn) && (index>=0));
res[i] += simuline[index];
}
}
} else { // TBM3, dim=3
Real dummy,halfz,ez;
assert(key->dim==3);
halfz = ((TBM_storage*)key->S)->half[2];
InvSqrtNlines=1.0 / sqrt((Real) TBM3D3_LINES);
for (n=0;n<TBM3D3_LINES;n++){
if ((GENERAL_PRINTLEVEL>=3) && (n % 100 == 0)) PRINTF("%d \n",n);
dummy = UNIFORM_RANDOM * PI;
ez = sin(dummy) * linescale;
dummy = cos(dummy) * linescale;
ey= UNIFORM_RANDOM * TWOPI;
ex= cos(ey) * dummy;
ey= sin(ey) * dummy;
simulate(key,simuline END_WITH_RANDOM);
offset= nnhalf - halfz * ez - halfy * ey - halfx * ex;
for (i=0;i<key->totallength;i++) {
register long index;
index =
(long) (offset + key->x[0][i] * ex + key->x[1][i] * ey + key->x[2][i] * ez);
assert((index<nn) && (index>=0));
res[i] += simuline[index];
}
}
}
}
}
{
register long i;
for(i=0;i<key->totallength;i++) { res[i] = res[i] * InvSqrtNlines; }
}
if (!(key->active=GENERAL_STORING)) {
assert(key->destruct!=NULL);
key->destruct(&key->S);
key->destruct=NULL;
}
return;
}