<< . .

. 48
( : 51)



. . >>

{violate=1;
if (height[iup[i]][j]>min)
{min=height[iup[i]][j];mini=iup[i];minj=j;}}
if (height[idown[i]][j]-height[i][j]>thresh)
{violate=1;
if (height[idown[i]][j]>min)
{min=height[idown[i]][j];mini=idown[i];minj=j;}}
if (height[i][jup[j]]-height[i][j]>thresh)
{violate=1;
if (height[i][jup[j]]>min)
{min=height[i][jup[j]];mini=i;minj=jup[j];}}
if (height[i][jdown[j]]-height[i][j]>thresh)
{violate=1;
if (height[i][jdown[j]]>min)
{min=height[i][jdown[j]];mini=i;minj=jdown[j];}}
if (violate==1)
{height[i][j]++;
height[mini][minj]--;
avalancheup(mini,minj);}
else
{hshadow=height[i][j]-0.5;
i2=iup[i];
while ((hshadow<mask[i2][j])&&(hshadow>0))
{if (height[i2][j]>=hshadow) mask[i2][j]=0; else mask[i2][j]=hshadow;
i2=iup[i2];hshadow-=0.5;}}
}

main()
{ FILE *fp1,*fp2;
int idum,t,l,ijump,jjump,i,j,duration;
float psand,pbed,p;
A6.1 WERNER™S EOLIAN DUNE MODEL 269


fp1=fopen("wernermodeltopo","w");
fp2=fopen("wernermodelshadowmask","w");
idum=-56;
psand=0.6;
pbed=0.4;
thresh=1;
l=5;
lattice_size_x=300;
lattice_size_y=300;
duration=100;
height=imatrix(1,lattice_size_x,1,lattice_size_y);
mask=matrix(1,lattice_size_x,1,lattice_size_y);
setupgridneighborsperiodic();
for (i=1;i<=lattice_size_x;i++)
for (j=1;j<=lattice_size_y;j++)
{height[i][j]=3;
mask[i][j]=0;}
for (t=1;t<=duration*lattice_size_x*lattice_size_y;t++)
{if (t%(lattice_size_x*lattice_size_y)==0) printf("%d\n",t);
ijump=(int)(ran3(&idum)*lattice_size_x)+1;
while (ijump>lattice_size_x) ijump=(int)(ran3(&idum)*lattice_size_x)+1;
jjump=(int)(ran3(&idum)*lattice_size_y)+1;
while (jjump>lattice_size_y) jjump=(int)(ran3(&idum)*lattice_size_y)+1;
while ((height[ijump][jjump]==0)||(mask[ijump][jjump]>0.1))
{ijump=(int)(ran3(&idum)*lattice_size_x)+1;
while (ijump>lattice_size_x) ijump=(int)(ran3(&idum)*lattice_size_x)+1;
jjump=(int)(ran3(&idum)*lattice_size_y)+1;
while (jjump>lattice_size_y)
jjump=(int)(ran3(&idum)*lattice_size_y)+1;}
height[ijump][jjump]--;
avalancheup(ijump,jjump);
ijump=(ijump+l)%lattice_size_x+1;
if (mask[ijump][jjump]>0.1) p=1;
else if (height[ijump][jjump]>0) p=psand; else p=pbed;
while (ran3(&idum)>p)
{ijump=(ijump+l)%lattice_size_x+1;
if (height[ijump][jjump]>0) p=psand; else p=pbed;}
height[ijump][jjump]++;
avalanchedown(ijump,jjump);}
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
{fprintf(fp1,"%d\n",height[i][j]);
fprintf(fp2,"%f\n",mask[i][j]);}
fclose(fp1);
fclose(fp2);
}
270 CODES FOR MODELING INSTABILITIES



A6.2 Oscillations in arid alluvial channels

The following code solves the equations for arid channel oscillations described in Section 7.5.

#define PI 3.141592653589793

main()
{ FILE *fp1;
int idum,*iup,*idown,count,printed,outputinterval,length_y,lattice_size_x,i;
float courantwave,courantdiff,maxdelw,del,*h2,*w2,h0,*flux,*deltaw,*deltah;
float w0,simulationlength,deltax,elapsedtime,timestep,channelslope;
float *height,*heightold,*width,*widthold,maxheightchange;

fp1=fopen("channeloscillations","w");
lattice_size_x=100;
deltax=100; /* m */
idum=-76;
channelslope=0.01; /* m/m */
w0=100; /* m */
h0=1; /* m */
iup=ivector(1,lattice_size_x);
idown=ivector(1,lattice_size_x);
height=vector(1,lattice_size_x);
heightold=vector(1,lattice_size_x);
deltah=vector(1,lattice_size_x);
width=vector(1,lattice_size_x);
widthold=vector(1,lattice_size_x);
deltaw=vector(1,lattice_size_x);
flux=vector(1,lattice_size_x);
h2=vector(1,lattice_size_x);
w2=vector(1,lattice_size_x);
simulationlength=5000000.0; /* yr */
outputinterval=100000.0;
for (i=1;i<=lattice_size_x;i++)
{width[i]=w0;
widthold[i]=width[i];
height[i]=-0.05*sin((i-1)*2*PI*3/(lattice_size_x-1))+0.03*gasdev(&idum);
heightold[i]=height[i];
flux[i]=0;h2[i]=0;w2[i]=0;
deltaw[i]=0;deltah[i]=0;
iup[i]=i+1;idown[i]=i-1;}
iup[lattice_size_x]=1;idown[1]=lattice_size_x;
timestep=1;
elapsedtime=0.0;
printed=0;
while (elapsedtime<simulationlength)
{printed=elapsedtime/outputinterval;
maxdelw=0;
A6.3 1D MODEL OF SPIRAL TROUGHS ON MARS 271


for (i=1;i<=lattice_size_x;i++)
flux[i]=((heightold[idown[i]]-heightold[i])/deltax+channelslope)/widthold[i];
for (i=1;i<=lattice_size_x;i++)
{del=0.5*timestep/deltax*(flux[idown[i]]-flux[i]);
h2[i]=0.5*(heightold[idown[i]]+heightold[i])+del;
w2[i]=widthold[i]-0.5*timestep*0.5*
(flux[idown[i]]+flux[i]-2*channelslope/w0)/(h0-h2[i]);}
for (i=1;i<=lattice_size_x;i++)
flux[i]=((h2[i]-h2[iup[i]])/deltax+channelslope)/w2[i];
for (i=1;i<=lattice_size_x;i++)
{del=timestep/deltax*(flux[i]-flux[iup[i]]);
height[i]+=del;
width[i]-=0.5*timestep*(flux[iup[i]]+flux[i]-
2*channelslope/w0)/(h0-heightold[i]);
if (fabs(widthold[idown[i]]-widthold[iup[i]])/widthold[i]>maxdelw)
maxdelw=fabs(widthold[idown[i]]-widthold[iup[i]])/widthold[i];
if (width[i]<10) width[i]=10;
if (width[i]>300) width[i]=300;}
elapsedtime+=timestep;
if ((maxdelw>0.5*deltax/timestep)||(timestep>deltax*deltax))
{elapsedtime-=timestep;
timestep/=2.0;
for (i=1;i<=lattice_size_x;i++)
{height[i]=heightold[i];
width[i]=widthold[i];}}
else
if ((maxdelw<0.5*deltax*deltax/timestep/10)&&(timestep<deltax*deltax/10))
timestep*=1.2;
for (i=1;i<=lattice_size_x;i++)
{heightold[i]=height[i];
widthold[i]=width[i];}
if ((int)(elapsedtime/outputinterval)>printed)
{for (i=1;i<=lattice_size_x;i++)
fprintf(fp1,"%d %f %f\n",i,height[i]/h0,width[i]/w0);
printed=(int)(elapsedtime/outputinterval);}}
fclose(fp1);
}


A6.3 1D model of spiral troughs on Mars

The following code integrates the 1D spiral trough model described in Section 7.7.

#define N 2
#define epssm 0.0001

float *xp,**yp,dxsav,a,eps;
int kmax,kount,*iup,*jup,*idown,*jdown;
272 CODES FOR MODELING INSTABILITIES


void derivs(x,y,dydx)
float dydx[],x,y[];
{
dydx[1]=y[2];
dydx[2]=(y[2]*(1-y[2])*(y[2]-a)-y[1])/eps;
}
main()
{ FILE *fp1;
int nbad,nok,t,lattice_size_x,i;
float *sol,h1=0.1,hmin=0.0,x1,x2,*lap,deltax,fac;
float *u,*v,timestep,kappa;



fp1=fopen("spiralmodel1d","w");
a=0.2;
eps=0.05;
deltax=0.2;
timestep=deltax*deltax/5;
kappa=timestep/(deltax*deltax);
lattice_size_x=512;
u=vector(1,lattice_size_x);
v=vector(1,lattice_size_x);
lap=vector(1,lattice_size_x);
sol=vector(1,2);
iup=ivector(1,lattice_size_x);
idown=ivector(1,lattice_size_x);
for (i=1;i<=lattice_size_x;i++)
{iup[i]=i+1;
idown[i]=i-1;}
iup[lattice_size_x]=1;
idown[1]=lattice_size_x;
for (i=1;i<=lattice_size_x;i++)
if ((i>lattice_size_x/2-5)&&(i<lattice_size_x/2-2))
{u[i]=0.25;v[i]=0.0;}
else
{u[i]=0;v[i]=0;}
/* plot initial condition first */
for (i=1;i<=lattice_size_x;i++)
fprintf(fp1,"%f %f %f\n",(i-lattice_size_x/2)/10.0,u[i],-v[i]);
fac=timestep*kappa/deltax/deltax;
x1=0;x2=timestep;
for (t=1;t<=1000;t++)
{for (i=1;i<=lattice_size_x;i++)
lap[i]=0;
for (i=1;i<=lattice_size_x;i++)
{sol[1]=v[i];sol[2]=u[i];
odeint(sol,N,x1,x2,epssm,h1,hmin,&nok,&nbad,derivs,rkqs);
v[i]=sol[1];u[i]=sol[2];
if (fabs(u[i])>0.00001) {lap[i]=-2*u[i];lap[iup[i]]=u[i];lap[idown[i]]=u[i];}}
A6.3 1D MODEL OF SPIRAL TROUGHS ON MARS 273


for (i=1;i<=lattice_size_x;i++)
if (fabs(lap[i])>0) u[i]+=fac*lap[i];}
for (i=1;i<=lattice_size_x;i++)
fprintf(fp1,"%f %f %f\n",(i-lattice_size_x/2)/10.0,u[i],-v[i]);
fclose(fp1);
}
Appendix 7




Codes for modeling stochastic processes

A7.1 Fractional-noise generation with Fourier-¬ltering method

The following code generates a 1D fractional noise (i.e. time series) using the Fourier-¬ltering method.
The user is prompted for all of the information used to generate the noise. The code will produce
noises with Gaussian or log-normal distributions depending on directions from the user. The random
number seed idum generates a pseudo-random sequence of numbers. In order to generate distinct
output, different values of idum must be used.


main()
{ float *precyear,beta,sum,sddevcomp,mean,sddev;
int length,year,format;
int idum;
char outputfile[30];
FILE *fp;

printf("\nName of output file: ");
scanf("%s",outputfile);
printf("\nLength of desired noise: (factor of two, please) ");
scanf("%d",&length);
printf("\nGaussian or Log-normal dataset? (1 for gaussian 2 for l-n): ");
scanf("%d",&format);
printf("\nMean: (if log-normal this mean is that of the final,
log transformed data) ");
scanf("%f",&mean);
printf("\nStandard Dev: ");
scanf("%f",&sddev);
printf("\nBeta: ");
scanf("%f",&beta);
printf("\nRandom number seed: (any negative integer) ");
scanf("%d",&idum);
fp=fopen(outputfile,"w");
precyear=vector(1,2*length);
/*we contruct a vector of length two times the input length so that we
may cut off the ends to eliminate the periodicity introduced by filtering*/
A7.1 FRACTIONAL-NOISE GENERATION WITH FOURIER-FILTERING METHOD 275


for (year=1;year<=2*length;year++)
precyear[year]=gasdev(&idum);
realft(precyear,length,1);
precyear[1]=0.0;
precyear[2]=precyear[2]/pow(0.5,beta/2.0);
for (year=3;year<=2*length;year++)
precyear[year]=precyear[year]/pow((year/2)/(float)(2*length),beta/2.0);
realft(precyear,length,-1);
sum=0.0;
for (year=length/2;year<=length+length/2;year++)
{sum+=precyear[year];}
sum=sum/length;
sddevcomp=0.0;
for (year=length/2;year<=length+length/2;year++)
sddevcomp+=(precyear[year]-sum)*(precyear[year]-sum);
sddevcomp=sddevcomp/length;
/*below we rescale the amplitudes to restore the desired moments*/
for (year=length/2;year<=length+length/2;year++)
{precyear[year]=(precyear[year]-sum)/sqrt(sddevcomp);
precyear[year]=(precyear[year]*sddev)+mean;
if (format==2) if (mean!=0.0) precyear[year]=
sqrt(log(1+sddev*sddev/(mean*mean)))*precyear[year]+

log(mean/sqrt(1+sddev*sddev/(mean*mean)));}
if (format==2) for (year=length/2;year<=length+length/2;year++)
precyear[year]=exp(precyear[year]);
for (year=length/2;year<=length+length/2;year++)
fprintf(fp,"%d %f\n",year-length/2,precyear[year]+(float)(year-length/2)/(length));
fclose(fp);
}

The following code generates a 2D Gaussian fractional noise (i.e. raster dataset) using the Fourier-
¬ltering method. The user is prompted for all of the information used to generate the dataset.
main()
{ float fact,*precyear,betax,sum,mean,sddev,sddevcomp;
int lengthx,lengthy;
int *nn,i,j,idum,index;
char outputfile[30];
FILE *fp;

printf("\nName of output file: ");
scanf("%s",outputfile);
printf("\nLength of desired noise: (factor of two, please) ");
scanf("%d",&lengthx);
lengthx=lengthx*2;
nn=ivector(1,2);
nn[1]=lengthx;
nn[2]=lengthx;
lengthy=lengthx;
276 CODES FOR MODELING STOCHASTIC PROCESSES


printf("\nMean: ");
scanf("%f",&mean);
printf("\nStandard Dev: ");
scanf("%f",&sddev);
printf("\nBeta: ");
scanf("%f",&betax);
printf("\nRandom number seed: (any negative integer) ");
scanf("%d",&idum);
fp=fopen(outputfile,"w");
precyear=vector(1,2*lengthx*lengthy);
for (i=1;i<=lengthx;i++)
for (j=1;j<=lengthy;j++)
{precyear[2*(i-1)*lengthy+2*j-1]=gasdev(&idum);
precyear[2*(i-1)*lengthy+2*j]=0.0;}
fourn(precyear,nn,2,1);
for (j=1;j<=lengthy;j++)
{precyear[2*j-1]=0.0;
precyear[2*j]=0.0;
precyear[2*(j-1)*lengthy+1]=0.0;
precyear[2*(j-1)*lengthy+2]=0.0;}
for (i=1;i<=lengthx/2;i++)
for (j=1;j<=lengthy/2;j++)
{fact=pow(sqrt(i*i+j*j)/(2*lengthx),(betax+1)/2.0);
index=2*i*lengthy+2*j+1;
precyear[index]=precyear[index]/(lengthx*lengthy*fact);
precyear[index+1]=precyear[index+1]/(lengthx*lengthy*fact);
index=2*i*lengthy+2*lengthy-2*j+1;
precyear[index]=precyear[index]/(lengthx*lengthy*fact);
precyear[index+1]=precyear[index+1]/(lengthx*lengthy*fact);
index=2*lengthx*lengthy-2*(i-1)*lengthy-2*(lengthy-j)+1;
precyear[index]=precyear[index]/(lengthx*lengthy*fact);
precyear[index+1]=precyear[index+1]/(lengthx*lengthy*fact);
index=2*lengthx*lengthy-2*lengthy-2*(i-1)*lengthy+2*(lengthy-j)+1;
precyear[index]=precyear[index]/(lengthx*lengthy*fact);
precyear[index+1]=precyear[index+1]/(lengthx*lengthy*fact);}
fourn(precyear,nn,2,-1);
sum=0.0;
for (i=1;i<=lengthx/2;i++)
for (j=1;j<=lengthy/2;j++)
sum+=precyear[2*(i-1)*lengthy+2*j-1];
sum=sum/(lengthx*lengthy);
sddevcomp=0.0;
for (i=1;i<=lengthx/2;i++)
for (j=1;j<=lengthy/2;j++)
sddevcomp+=(precyear[2*(i-1)*lengthy+2*j-1]-sum)*
(precyear[2*(i-1)*lengthy+2*j-1]-sum);
sddevcomp=sqrt(sddevcomp/(lengthx*lengthy));
for (i=1;i<=lengthx/2;i++)
for (j=1;j<=lengthy/2;j++)
A7.2 STOCHASTIC MODEL OF PLEISTOCENE ICE AGES 277


{index=2*(i-1)*lengthy+2*j-1;
precyear[index]=(precyear[index]-mean)/sddevcomp;
precyear[index]=(precyear[index]*sddev)+mean;
fprintf(fp,"%f\n",precyear[index]);}
fclose(fp);
}


A7.2 Stochastic model of Pleistocene ice ages

The following code is the basis for the stochastic-resonance model of ice ages. The model outputs
four columns with the form: t (yr), T , v , and B h .
#define maxoutputlength 10000000

main()
{ FILE *fp1;

<< . .

. 48
( : 51)



. . >>