/*
    HWT_3D_mm: Acoustical 3-D Huygens wavefront/ray tracing for moving media.
    An extension for moving media by N. Zabotin, O.A. Godin and L.Ye. Zabotina
    (Copyright (C) 2013 University of Colorado Boulder) of an algorithm
    developed by P. Sava and S. Fomel as a part of Madagascar open project
    (Copyright (C) 2007 Colorado School of Mines).

    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 (file COPYING.txt).

    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

Make file for Linux (tested with Intel C++ compiler): HWT_3D_mm_make
For Windows: Use MS Visual Studio standard procedures (tested with Intel compiler)

USAGE:
     HWT_3D_mm par=<name of text file containing list of parameters
                     for specific mode of calculation; example: par_timefront.txt>
PARAMETERS:
infile - name of text file containing parameters of the model of propagation medium:
          description of the grid used and names of binary files with sound speed and
          current data; example files for a model with linear profiles: linz3d_u1.txt,
          a01_u05_gama005.rsf, a01_u05_gama005_x.rsf, a01_u05_gama005_y.rsf; example
          program write_rsf_3d.cpp demonstrates creation of the binary files.
dt,nt - time step (in s) and total number of time steps;
ng,og,dg - number of steps, initial value, and step for the grazing angle (in deg);
nh,oh,dh - number of steps, initial value, and step for the bearing angle (in deg);
xsou,ysou,zsou - coordinates of a point source (in km);
refl - switch reflections from the plane surface z=0 on (=y) or off (=n).

Input parameters triggering specific modes of calculation:
1. 1D timefront calculation (time dependence of 1D coordinates of intersection
                          of a 3D wavefront with a straight line set in 3D
                          space by coordinates of two points on it)
tR - a set of six floating numbers (x,y,z coordinates of the 1st point and
                                    x,y,z coordinates of the 2nd point,
                                    separated by commas (in km)).

2. Wavefront calculation (3D shape for a specific time)
tw - specific time for wavefront calculation (in s).

3. Intersection of a 3D wavefront with horizontal (z=0) plane at a specific time
tz - specific time for the wavefront intersection calculation (in s).

4. Intersection of a 3D wavefront with vertical (y=0) plane at a specific time
ty - specific time for the wavefront intersection calculation (in s).

5. 3D ray trajectories x(t),y(t),z(t) are written in output file if no parameters
mentioned in items 1-4 are set.

*/

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <sys/types.h>

#include "alloc.h"
#include "_bool.h"
#include "axa.h"
#include "_defs.h"
#include "error.h"
#include "file.h"
#include "getpar.h"
#include "point.h"
#include "simtab.h"
#include "vector.h"
#include "HWT_3D_mm.h"

int main(int argc, char* argv[])
{
    bool verb,tfront=false,refl=false;
    float        scaleray;
    sf_axis az,ax,ay;    /* Cartesian coordinates */
    sf_axis at,ag,ah,aa; /* Ray coordinates */
    int it,ig,ih;

    float xsou,ysou,zsou; /* source coordinates */

    sf_file Fv; /*   velocity file */

    float ***vv=NULL; /* velocity       */
    float ***curx=NULL; /* current x-component*/
    float ***cury=NULL; /* current z-component*/
    float ***curz=NULL; /* current z-component*/

    pt3d_cur   **wm=NULL; /* wavefront it-1 */
    pt3d_cur   **wo=NULL; /* wavefront it   */
    pt3d_cur   **wp=NULL; /* wavefront it+1 */
    pt3d_cur Tm,To,Tp; /* points on ray ig,ih:    it-1,it,it+1 */

    int n;
    float o,d;
    float tt=0;
    vm   v;
    int smth=1;
    char* tag;

    double tw,tz,ty;
// z reflections
    pt3d_cur   **wg=NULL;
    bool       **lg=NULL;
    float Asurf=0.0, Bsurf=0.0, Csurf=1.0, dsurf=0.0; /* horizontal plane*/
    point3d pz0;

// y reflections
    double Apy=0.0, Bpy=1.0, Cpy=0.0, Dpy=0;  /* vertical plane*/
    point3d py0;

// timefront
    ray3d R={0.0};
    triangle3d T;
    point3d Intsec;
    int cod=0;
    double dd,g,h;
    size_t tn=6;
    double *tR;

    FILE *stream;
    char*outfile;

/*------------------------------------------------------------*/

    sf_init(argc,argv);

/*Parameters from command line*/
    if(! sf_getbool("verb",&verb))     verb=false; // output on the screen y/n
    if(! sf_getfloat("scaleray",&scaleray)) scaleray=1.0;
    if(! sf_getbool("refl",&refl))     refl=false; // surface reflections y/n

    tR=sf_doublealloc(tn);
    if ( sf_getdoubles("tR",tR,tn)) {
      tfront=true;
      memcpy(&R,&tR[0],6*sizeof(double));
     }
    free(tR);

    if(! sf_getdouble("tw",&tw))   tw=0.0;   // Wavefront calculation
    if(! sf_getdouble("ty",&ty))   ty=0.0;   // Intersection of a 3D wavefront with vertical (y=0) plane
    if(! sf_getdouble("tz",&tz))   tz=0.0;   // Intersection of a 3D wavefront with horizontal (z=0) plane

    /* input file */
    tag=sf_getstring("infile");
    Fv = sf_input( tag );  // ("in");

    /* information about axes from input file */
    az=sf_iaxa(Fv,1); sf_setlabel(az,"z");
    ax=sf_iaxa(Fv,2); sf_setlabel(ax,"x");
    ay=sf_iaxa(Fv,3); sf_setlabel(ay,"y");

   if(verb) {  /*print axes information to the screen*/
     sf_raxa(az);
     sf_raxa(ax);
     sf_raxa(ay);

    fprintf(stderr,"tw=%f tfront=%d ty=%f tz=%f\n",tw,tfront,ty,tz);
    fprintf(stderr,"R=%f %f %f %f %f %f\n",R.P0.x,R.P0.y,R.P0.z,R.P1.x,R.P1.y,R.P1.z);
    fprintf(stderr,"refl=%d\n",refl);
    }

    /* velocity file */
    vv=sf_floatalloc3(sf_n(az),sf_n(ax),sf_n(ay));
    sf_floatread(vv[0][0],sf_n(az)*sf_n(ax)*sf_n(ay),Fv);
    if (verb) fprintf(stderr,"vv0=%f vv1=%f \n",vv[0][0][0],vv[6][10][99]);

    /* current file */
    curx=sf_floatalloc3(sf_n(az),sf_n(ax),sf_n(ay));
    cury=sf_floatalloc3(sf_n(az),sf_n(ax),sf_n(ay));
    curz=sf_floatalloc3(sf_n(az),sf_n(ax),sf_n(ay));
	for( it=0; it<sf_n(ay); it++) {
	for( ig=0; ig<sf_n(ax); ig++) {
    for( ih=0; ih<sf_n(az); ih++) {
		curz[it][ig][ih]=curx[it][ig][ih]=cury[it][ig][ih]=0.0;
   	}}}
    cod=current_read (curx[0][0],cury[0][0],curz[0][0],sf_n(az)*sf_n(ax)*sf_n(ay),Fv);
    if (cod==-1)    exit(0);
    else
	if (verb) fprintf(stderr,"curx=%f cury=%f curz=%f \n",curx[0][0][0],cury[0][0][0],curz[0][0][0]);


    /* source location */
    if(! sf_getfloat("xsou",&xsou)) xsou=sf_o(ax) + sf_n(ax)*sf_d(ax)/2;
    if(! sf_getfloat("ysou",&ysou)) ysou=sf_o(ay) + sf_n(ay)*sf_d(ay)/2;
    if(! sf_getfloat("zsou",&zsou)) zsou=sf_o(az) + sf_n(az)*sf_d(az)/2;
    if(verb) sf_warning("xsou=%f ysou=%f zsou=%f",xsou,ysou,zsou);

    /* time axis */
    if(! sf_getint  ("nt",&n)) n=100;
    if(! sf_getfloat("ot",&o)) o=0;
    if(! sf_getfloat("dt",&d)) d=0.001;
    at=sf_maxa(n,o,d);
    if (verb) fprintf(stderr,"ot=%f dt=%f at=%d\n",sf_o(at),sf_d(at),sf_n(at));

    /* shooting angle axis */
    if(! sf_getint  ("ng",&n)) n=360;
    if(! sf_getfloat("og",&o)) o=-180;
    if(! sf_getfloat("dg",&d)) d=1;
    ag=sf_maxa(n,o,d);

    if(! sf_getint  ("nh",&n)) n=360;
    if(! sf_getfloat("oh",&o)) o=-180;
    if(! sf_getfloat("dh",&d)) d=1;
    ah=sf_maxa(n,o,d);
    if(sf_n(ay)==1) {  sf_setn(ah,1); sf_seto(ah,0); } /* 2-D case */

	if (verb) {
    fprintf(stderr,"ag %f %f %d\n",sf_o(ag),sf_d(ag),sf_n(ag));
    fprintf(stderr,"ah %f %f %d\n",sf_o(ah),sf_d(ah),sf_n(ah));
	}
    /*------------------------------------------------------------*/
    /* Open text file to write wavefronts or timefronts*/
    if (tfront)  outfile="timefront.txt";
    else         outfile="wavefront.txt";
    if ( (stream  = fopen( outfile, "w+t" )) == NULL  )  fprintf(stderr,"%s\n","Can not open output file!");
    /*------------------------------------------------------------*/

    /* allocate wavefronts */
    wm = pt3d_curalloc2(sf_n(ag),sf_n(ah));
    wo = pt3d_curalloc2(sf_n(ag),sf_n(ah));
    wp = pt3d_curalloc2(sf_n(ag),sf_n(ah));

    wg = pt3d_curalloc2(sf_n(ag),sf_n(ah));
    lg = sf_boolalloc2(sf_n(ag),sf_n(ah));


    /* initialize wavefronts */
    for( ih=0; ih<sf_n(ah); ih++) {
	for( ig=0; ig<sf_n(ag); ig++) {
	    wm[ih][ig].x=wo[ih][ig].x=wp[ih][ig].x=wg[ih][ig].x=0;
	    wm[ih][ig].y=wo[ih][ig].y=wp[ih][ig].y=wg[ih][ig].y=0;
	    wm[ih][ig].z=wo[ih][ig].z=wp[ih][ig].z=wg[ih][ig].z=0;
	    wm[ih][ig].v=wo[ih][ig].v=wp[ih][ig].v=wg[ih][ig].v=0;

        wm[ih][ig].mx = wo[ih][ig].mx = wp[ih][ig].mx=wg[ih][ig].mx=0;
        wm[ih][ig].my = wo[ih][ig].my = wp[ih][ig].my=wg[ih][ig].my=0;
        wm[ih][ig].mz = wo[ih][ig].mz = wp[ih][ig].mz=wg[ih][ig].mz=0;

        lg[ih][ig]=false;
	}
    }

/*------------------------------------------------------------*/
    /* init HWT */
     hwt3d_init(vv,curx,cury,curz,az,ax,ay,at,ag,ah,Asurf,Bsurf,Csurf,dsurf);
/*------------------------------------------------------------*/

    /* construct it=0 wavefront */
    it=0;
    if(verb) sf_warning("it=%d of %d",it+1,sf_n(at));
    for( ih=0; ih<sf_n(ah); ih++) {
	for( ig=0; ig<sf_n(ag); ig++) {
	    wm[ih][ig].x=xsou;
	    wm[ih][ig].y=ysou;
	    wm[ih][ig].z=zsou;
	    v=hwt3d_getv(wm[ih][ig]);
        wm[ih][ig]=get_point(v,wm[ih][ig]);
	}
    }
    if (!tw && !tfront && !tz && !ty){
       tt=sf_o(at)+it*sf_d(at);
       wavefront(stream,sf_n(ah),sf_n(ag),wm, sf_o(ah),sf_o(ag),sf_d(ah),sf_d(ag),tt);
    }
/*------------------------------------------------------------*/

    /* construct it=1 wavefront */
    it=1;
    if(verb) sf_warning("it=%d of %d",it+1,sf_n(at));
    for( ih=0; ih<sf_n(ah); ih++) {
	for( ig=0; ig<sf_n(ag); ig++) {

	    dd = sf_d(at);
	    g = (sf_o(ag)+ig*sf_d(ag)) * SF_PI/180;
	    h = (sf_o(ah)+ih*sf_d(ah)) * SF_PI/180;

	    wo[ih][ig].x=xsou + dd*sin(g)*cos(h)*wm[ih][ig].v + dd*wm[ih][ig].mx;
	    wo[ih][ig].y=ysou + dd*sin(g)*sin(h)*wm[ih][ig].v + dd*wm[ih][ig].my;
	    wo[ih][ig].z=zsou + dd*cos(g)*wm[ih][ig].v        + dd*wm[ih][ig].mz;
	    v=hwt3d_getv(wo[ih][ig]);
        wo[ih][ig]=get_point(v,wo[ih][ig]);
	}
    }
    if (!tw && !tfront && !tz && !ty){
       tt=sf_o(at)+it*sf_d(at);
       wavefront(stream,sf_n(ah),sf_n(ag),wo, sf_o(ah),sf_o(ag),sf_d(ah),sf_d(ag),tt);
    }

/*------------------------------------------------------------*/
    /* LOOP over time */
    for (it=2; it<sf_n(at); it++) {
        smth++;
        tt=sf_o(at)+it*sf_d(at);
        sf_warning("it=%d of %d, tt=%f",it+1,sf_n(at),tt);

	for (ih=0; ih<sf_n(ah); ih++) {
    for (ig=0; ig<sf_n(ag); ig++) {

        To = wo[ih][ig];    // t wavefront
		if (lg[ih][ig]) {
                   Tm = wg[ih][ig];
                   lg[ih][ig]=false;
        }
        else Tm = wm[ih][ig];  // t-1 wavefront
        Tp = hwt3d_raytr(Tm,To,scaleray,stream,refl);


//------------------------   Reflections   ---------------------------------
    if ( lg_surf && refl) {  // from horizontal surface

       double pp=(Asurf*To.x+Bsurf*To.y+Csurf*To.z+dsurf)/(Asurf*Asurf+Bsurf*Bsurf+Csurf*Csurf);
       wg[ih][ig].x=To.x-2.0*Asurf*pp;
       wg[ih][ig].y=To.y-2.0*Bsurf*pp;
       wg[ih][ig].z=To.z-2.0*Csurf*pp;
       v=hwt3d_getv(wg[ih][ig]);   wg[ih][ig]=get_point(v,wg[ih][ig]);
       lg[ih][ig]=true;

//looking for intersection of the wavefront with z=0
// So - the point should have been reflected,   Tp - new reflected point
    if ( tz && tt>=tz )  {
       pz0=zeropoint(0);
       pz0 = zinterp(To,So);  //Tp
       g=sf_o(ag)+ig*sf_d(ag);
       h=sf_o(ah)+ih*sf_d(ah);
       fprintf(stream,"%8.2f %8.2f %10.5f %10.5f %6.1f\n",g,h,pz0.x,pz0.y,tt);
       So=zeropoint(0);
    }
    }

//--------------------------------------------------------------------------
		wp[ih][ig] = Tp;
	    }
	}

	/* write wavefront */
        if (!ty && !tfront && !tz)
        if (tw) {
        if (fabs(tt-tw)<sf_d(at)/10) wavefront(stream,sf_n(ah),sf_n(ag),wp, sf_o(ah),sf_o(ag),sf_d(ah),sf_d(ag),tt);
        }
        else wavefront(stream,sf_n(ah),sf_n(ag),wp, sf_o(ah),sf_o(ag),sf_d(ah),sf_d(ag),tt);



//---------------------  y-plane intersections ---------------------
    if (ty && !tfront && !tz) {
    for (ih=0; ih<sf_n(ah)-1; ih++) {
    for (ig=0; ig<sf_n(ag)-1; ig++) {

        if ( (Apy*wp[ih][ig].x+Bpy*wp[ih][ig].y+Cpy*wp[ih][ig].z+Dpy)*(Apy*wp[ih][ig+1].x+Bpy*wp[ih][ig+1].y+Cpy*wp[ih][ig+1].z+Dpy)<0 ){
           if ( fabs(tt-ty)<sf_d(at)/10 )  {
              py0=zeropoint(0);
              py0 = yinterp(wp[ih][ig],wp[ih][ig+1],Dpy);
              g=sf_o(ag)+ig*sf_d(ag);         h=sf_o(ah)+ih*sf_d(ah);
              fprintf(stream,"%8.2f %8.2f %10.5f %10.5f %10.5f %6.1f\n",g,h,py0.x,py0.y,py0.z,tt);
           }
        }
        if ( (Apy*wp[ih][ig].x+Bpy*wp[ih][ig].y+Cpy*wp[ih][ig].z+Dpy)*(Apy*wp[ih+1][ig].x+Bpy*wp[ih+1][ig].y+Cpy*wp[ih+1][ig].z+Dpy)<0 ){
           if ( fabs(tt-ty)<sf_d(at)/10 )  {
              py0=zeropoint(0);
              py0 = yinterp(wp[ih][ig],wp[ih+1][ig],Dpy);
              g=sf_o(ag)+ig*sf_d(ag);        h=sf_o(ah)+ih*sf_d(ah);
              fprintf(stream,"%8.2f %8.2f %10.5f %10.5f %10.5f %6.1f\n",g,h,py0.x,py0.y,py0.z,tt);
           }
        }

    }}
    }

	/* step in time */
	for( ih=0; ih<sf_n(ah); ih++) {
	    for( ig=0; ig<sf_n(ag); ig++) {
		wm[ih][ig] = wo[ih][ig];
		wo[ih][ig] = wp[ih][ig];

        /* timefront */
        if (tfront && !ty && !tw && !tz) {
                if (ih>0 && ig<sf_n(ag)-1 && tt>=5) {
                   T.V0=get_point3d(wp[ih-1][ig]);
                   T.V1=get_point3d(wp[ih][ig]);
                   T.V2=get_point3d(wp[ih][ig+1]);

                   cod=intersect_RayTriangle( R, T, &Intsec);
                   g=sf_o(ag)+ig*sf_d(ag);
                   h=sf_o(ah)+ih*sf_d(ah);
                   if (cod==1) fprintf(stream,"%7.3f %7.3f %10.5f %10.5f %10.5f %10.5f\n",g,h,tt,Intsec.x,Intsec.y,Intsec.z);
                }
        }

	    }
	}
    } /* end it */
    fclose( stream );
    /*------------------------------------------------------------*/

    free(**vv); free(*vv); free(vv);
    ;           free(*wm); free(wm);
    ;           free(*wo); free(wo);
    ;           free(*wp); free(wp);

    /*------------------------------------------------------------*/


    exit (0);
}