diff --git a/fdelmodc/decomposition.c.new b/fdelmodc/decomposition.c.new
deleted file mode 100644
index 3b03ba5a6c5956578ee28c89977bdde47c9d073f..0000000000000000000000000000000000000000
--- a/fdelmodc/decomposition.c.new
+++ /dev/null
@@ -1,247 +0,0 @@
-/*
- * decomposition.c
- *
- * Kees Wapenaar "Reciprocity properties of one-way propagators"
- * GEOPHYSICS, VOL. 63, NO. 4 (JULY-AUGUST 1998); P. 1795–1798
- *
- * Created by Jan Thorbecke on 17/03/2014.
- * Copyright 2014 TU Delft. All rights reserved.
- *
- */
-
-#include <assert.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <errno.h>
-#include <math.h>
-#include <string.h>
-
-#define MAX(x,y) ((x) > (y) ? (x) : (y))
-#define MIN(x,y) ((x) < (y) ? (x) : (y))
-
-#ifndef COMPLEX
-typedef struct _complexStruct { /* complex number */
- float r,i;
-} complex;
-typedef struct _dcomplexStruct { /* complex number */
- double r,i;
-} dcomplex;
-#endif/* complex */
-
-complex firoot(float x, float stab)
-complex ciroot(complex x, float stab);
-complex cwp_csqrt(complex z);
-
-void decud(float om, float rho, float cp, float dx, int nkx, complex *pu);
-
-void kxwfilter(complex *data, float k, float dx, int nkx,
- float alfa1, float alfa2, float perc);
-
-void kxwdecomp(complex *rp, complex *rvz, complex *up, complex *down,
- int nkx, float dx, int nt, float dt, float fmin, float fmax,
- float cp, float rho)
-{
- int iom, iomin, iomax, ikx, nfreq;
- float omin, omax, deltom, om, kp, df;
- float alpha, eps;
- complex *pu, w;
- complex ax, az;
-
- df = 1.0/((float)nt*dt);
- deltom = 2.*M_PI*df;
- omin = 2.*M_PI*fmin;
- omax = 2.*M_PI*fmax;
- nfreq = nt/2+1;
- eps = 0.01;
- alpha = 0.1;
-
- iomin = (int)MIN((omin/deltom), (nfreq-1));
- iomin = MAX(iomin, 1);
- iomax = MIN((int)(omax/deltom), (nfreq-1));
-
- pu = (complex *)malloc(nkx*sizeof(complex));
-
- for (iom = iomin; iom <= iomax; iom++) {
- om = iom*deltom;
-
- decud(om, rho, cp, dx, nkx, alpha, eps, pu);
-
-/*
- kxwfilter(dpux, kp, dx, nkx, alfa1, alfa2, perc);
- kxwfilter(dpuz, kp, dx, nkx, alfa1, alfa2, perc);
-*/
-
- for (ikx = 0; ikx < nkx; ikx++) {
- ax.r = 0.5*rp[iom*nkx+ikx].r;
- ax.i = 0.5*rp[iom*nkx+ikx].i;
- az.r = 0.5*(rvz[iom*nkx+ikx].r*pu[ikx].r - rvz[iom*nkx+ikx].i*pu[ikx].i);
- az.i = 0.5*(rvz[iom*nkx+ikx].i*pu[ikx].r + rvz[iom*nkx+ikx].r*pu[ikx].i);
-
- down[iom*nkx+ikx].r = ax.r + az.r;
- down[iom*nkx+ikx].i = ax.i + az.i;
- up[iom*nkx+ikx].r = ax.r - az.r;
- up[iom*nkx+ikx].i = ax.i - az.i;
- }
-
- }
-
- free(pu);
-
- return;
-}
-
-void decud(float om, float rho, float cp, float dx, int nkx, float alpha, float eps, complex *pu);
-{
- int ikx, ikxmax1, ikxmax2, filterpoints, filterppos;
- float mu, kp, kp2, ks, ks2, ksk;
- float kx, kx2, kzp2, kzs2, dkx;
- float kxfmax, kxnyq, kpos, kneg, alfa, kfilt, perc, band, *filter;
- complex kzp, kzs, cste, ckp, ckp2, ckzp2;
-
-/* with complex frequency
- wom.r=om;
- wom.i=alpha;
-
- ckp.r = wom.r/cp;
- ckp.i = wom.i/cp;
- ckp2.r = ckp.r*ckp.r-ckp.i*ckp.i;
- ckp2.i = 2.0*ckp.r*ckp.i;
- stab = eps*eps*(ckp.r*ckp.r+ckp.i*ckp.i);
-*/
-
-
- kp = om/cp;
- kp2 = kp*kp;
- dkx = 2.0*M_PI/(nkx*dx);
- stab = eps*eps*kp*kp;
-
- /* make kw filter at maximum angle alfa */
- alfa = 90.0;
- perc = 0.10; /* percentage of band to use for smooth filter */
- filter = (float *)malloc(nkx*sizeof(float));
- kpos = kp*sin(M_PI*alfa/180.0);
- kneg = -kpos;
- kxnyq = M_PI/dx;
- if (kpos > kxnyq) kpos = kxnyq;
- band = kpos;
- filterpoints = (int)fabs((int)(perc*band/dkx));
- kfilt = fabs(dkx*filterpoints);
- if (kpos+kfilt < kxnyq) {
- kxfmax = kpos+kfilt;
- filterppos = filterpoints;
- }
- else {
- kxfmax = kxnyq;
- filterppos = (int)(0.15*nkx/2);
- }
- ikxmax1 = (int) (kxfmax/dkx);
- ikxmax2 = ikxmax1 - filterppos;
- // fprintf(stderr,"ikxmax1=%d ikxmax2=%d nkp=%d nkx=%d\n", ikxmax1, ikxmax2, (int)(kp/dkx), nkx);
-
- for (ikx = 0; ikx < ikxmax2; ikx++)
- filter[ikx]=1.0;
- for (ikx = ikxmax2; ikx < ikxmax1; ikx++)
- filter[ikx] =(cos(M_PI*(ikx-ikxmax2)/(ikxmax1-ikxmax2))+1)/2.0;
- for (ikx = ikxmax1; ikx <= nkx/2; ikx++)
- filter[ikx] = 0.0;
- /* end of kxfilter */
-
- for (ikx = 0; ikx <= (nkx/2); ikx++) {
- kx = ikx*dkx;
- kx2 = kx*kx;
- kzp2 = kp2 - kx2;
- kzp = firoot(kzp2);
-
-/* with complex frequency
- kzp2.r = kp2.r - kx2;
- kzp2.i = kp2.i;
- kzp = ciroot(kzp2, stab);
-*/
- if (kzp2 != 0) {
- pu[ikx].r = filter[ikx]*om*rho*kzp.r;
- pu[ikx].i = filter[ikx]*om*rho*kzp.i;
-// pu[ikx].r = om*rho*kzp.r;
-// pu[ikx].i = om*rho*kzp.i;
- }
- else {
- pu[ikx].r = 0.0;
- pu[ikx].i = 0.0;
- }
- }
-
-/* operators are symmetric in kx-w domain */
- for (ikx = (nkx/2+1); ikx < nkx; ikx++) {
- pu[ikx] = pu[nkx-ikx];
- }
- free(filter);
-
- return;
-}
-
-/* compute 1/x */
-complex firoot(float x, float stab)
-{
- complex z;
- if (x > 0.0) {
- z.r = 1.0/sqrt(x+stab);
- z.i = 0.0;
- }
- else if (x == 0.0) {
- z.r = 0.0;
- z.i = 0.0;
- }
- else {
- z.r = 0.0;
- z.i = 1.0/sqrt(-x+stab);
- }
- return z;
-}
-
-complex ciroot(complex x, float stab)
-{
- complex z, kz, kzz;
- float kd;
-
- if (x.r == 0.0) {
- z.r = 0.0;
- z.i = 0.0;
- }
- else {
- kzz = cwp_csqrt(x);
- kz.r = kzz.r;
- kz.i = -abs(kzz.i);
- kd = kz.r*kz.r+kz.i*kz.i+stab;
- z.r = kz.r/kd;
- z.i = -kz.i/kd;
- }
- return z;
-}
-
-complex cwp_csqrt(complex z)
-{
- complex c;
- float x,y,w,r;
- if (z.r==0.0 && z.i==0.0) {
- c.r = c.i = 0.0;
- return c;
- } else {
- x = fabs(z.r);
- y = fabs(z.i);
- if (x>=y) {
- r = y/x;
- w = sqrt(x)*sqrt(0.5*(1.0+sqrt(1.0+r*r)));
- } else {
- r = x/y;
- w = sqrt(y)*sqrt(0.5*(r+sqrt(1.0+r*r)));
- }
- if (z.r>=0.0) {
- c.r = w;
- c.i = z.i/(2.0*w);
- } else {
- c.i = (z.i>=0.0) ? w : -w;
- c.r = z.i/(2.0*c.i);
- }
- return c;
- }
-}
-
diff --git a/fdelmodc/writeSnapTimes.c.org b/fdelmodc/writeSnapTimes.c.org
deleted file mode 100644
index 69120ea518a619ede9580aec09e0c92f8f42a26c..0000000000000000000000000000000000000000
--- a/fdelmodc/writeSnapTimes.c.org
+++ /dev/null
@@ -1,172 +0,0 @@
-#define _FILE_OFFSET_BITS 64
-#define _LARGEFILE_SOURCE
-#define _LARGEFILE64_SOURCE
-
-#include <assert.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <errno.h>
-#include <math.h>
-#include <string.h>
-#include "par.h"
-#include "segy.h"
-#include "fdelmodc.h"
-
-/**
-* Writes gridded wavefield(s) at a desired time to output file(s)
-*
-* AUTHOR:
-* Jan Thorbecke (janth@xs4all.nl)
-* The Netherlands
-**/
-
-
-FILE *fileOpen(char *file, char *ext, int append);
-int traceWrite(segy *hdr, float *data, int n, FILE *fp);
-
-#define MAX(x,y) ((x) > (y) ? (x) : (y))
-#define MIN(x,y) ((x) < (y) ? (x) : (y))
-#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
-
-int writeSnapTimes(modPar mod, snaPar sna, int ixsrc, int izsrc, int itime, float *vx, float *vz, float *tzz, float *txx, float *txz, int verbose)
-{
- FILE *fpvx, *fpvz, *fptxx, *fptzz, *fptxz, *fpp, *fppp, *fpss;
- int append, isnap;
- int n1, ibnd, ixs, izs, ize, i, j;
- int ix, iz, ix2, iz2;
- float *snap, sdx;
- segy hdr;
-
- if (sna.nsnap==0) return 0;
-
- ibnd = mod.iorder/2-1;
- n1 = mod.naz;
- sdx = 1.0/mod.dx;
-
- /* check if this itime is a desired snapshot time */
- if ( (((itime-sna.delay) % sna.skipdt)==0) &&
- (itime >= sna.delay) &&
- (itime <= sna.delay+(sna.nsnap-1)*sna.skipdt) ) {
-
- isnap = NINT((itime-sna.delay)/sna.skipdt);
- if (verbose) vmess("Writing snapshot(%d) at time=%.3f", isnap+1, itime*mod.dt);
-
- if (isnap) append=1;
- else append=0;
-
- if (sna.type.vx) fpvx = fileOpen(sna.file_snap, "_svx", append);
- if (sna.type.vz) fpvz = fileOpen(sna.file_snap, "_svz", append);
- if (sna.type.p) fpp = fileOpen(sna.file_snap, "_sp", append);
- if (sna.type.txx) fptxx = fileOpen(sna.file_snap, "_stxx", append);
- if (sna.type.tzz) fptzz = fileOpen(sna.file_snap, "_stzz", append);
- if (sna.type.txz) fptxz = fileOpen(sna.file_snap, "_stxz", append);
- if (sna.type.pp) fppp = fileOpen(sna.file_snap, "_spp", append);
- if (sna.type.ss) fpss = fileOpen(sna.file_snap, "_sss", append);
-
- memset(&hdr,0,TRCBYTES);
- hdr.dt = 1000000*(mod.dt);
- hdr.scalco = -1000;
- hdr.scalel = -1000;
- hdr.sx = 1000*(mod.x0+ixsrc*mod.dx);
- hdr.sdepth = 1000*(mod.z0+izsrc*mod.dz);
- hdr.fldr = isnap+1;
- hdr.trid = 1;
- hdr.ns = sna.nz;
- hdr.trwf = sna.nx;
- hdr.ntr = (isnap+1)*sna.nx;
- hdr.f1 = sna.z1*mod.dz+mod.z0;
- hdr.f2 = sna.x1*mod.dx+mod.x0;
- hdr.d1 = mod.dz*sna.skipdz;
- hdr.d2 = mod.dx*sna.skipdx;
-
-/***********************************************************************
-* vx velocities have one sample less in x-direction
-* vz velocities have one sample less in z-direction
-* txz stresses have one sample less in z-direction and x-direction
-***********************************************************************/
-
- snap = (float *)malloc(sna.nz*sizeof(float));
-
- /* Decimate, with skipdx and skipdz, the number of gridpoints written to file
- and write to file. */
- for (ixs=sna.x1, i=0; ixs<=sna.x2; ixs+=sna.skipdx, i++) {
- hdr.tracf = i+1;
- hdr.tracl = isnap*sna.nx+i+1;
- hdr.gx = 1000*(mod.x0+ixs*mod.dx);
- ix = ixs+ibnd;
- ix2 = ix+1;
-
- izs = sna.z1+ibnd;
- ize = sna.z2+ibnd;
- if (sna.type.vx) {
- for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
- snap[j] = vx[ix2*n1+iz];
- }
- traceWrite(&hdr, snap, sna.nz, fpvx);
- }
- if (sna.type.vz) {
- for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
- snap[j] = vz[ix*n1+iz+1];
- }
- traceWrite(&hdr, snap, sna.nz, fpvz);
- }
- if (sna.type.p) {
- for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
- snap[j] = tzz[ix*n1+iz];
- }
- traceWrite(&hdr, snap, sna.nz, fpp);
- }
- if (sna.type.tzz) {
- for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
- snap[j] = tzz[ix*n1+iz];
- }
- traceWrite(&hdr, snap, sna.nz, fptzz);
- }
- if (sna.type.txx) {
- for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
- snap[j] = txx[ix*n1+iz];
- }
- traceWrite(&hdr, snap, sna.nz, fptxx);
- }
- if (sna.type.txz) {
- for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
- snap[j] = txz[ix2*n1+iz+1];
- }
- traceWrite(&hdr, snap, sna.nz, fptxz);
- }
- /* calculate divergence of velocity field */
- if (sna.type.pp) {
- for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
- iz2 = iz+1;
- snap[j] = sdx*((vx[ix2*n1+iz]-vx[ix*n1+iz])+
- (vz[ix*n1+iz2]-vz[ix*n1+iz]));
- }
- traceWrite(&hdr, snap, sna.nz, fppp);
- }
- /* calculate rotation of velocity field */
- if (sna.type.ss) {
- for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
- iz2 = iz+1;
- snap[j] = sdx*((vx[ix2*n1+iz2]-vx[ix2*n1+iz])-
- (vz[ix2*n1+iz2]-vz[ix*n1+iz2]));
- }
- traceWrite(&hdr, snap, sna.nz, fpss);
- }
-
- }
-
- if (sna.type.vx) fclose(fpvx);
- if (sna.type.vz) fclose(fpvz);
- if (sna.type.p) fclose(fpp);
- if (sna.type.txx) fclose(fptxx);
- if (sna.type.tzz) fclose(fptzz);
- if (sna.type.txz) fclose(fptxz);
- if (sna.type.pp) fclose(fppp);
- if (sna.type.ss) fclose(fpss);
-
- free(snap);
- }
-
- return 0;
-}
-