#include <stdio.h>
#include <math.h>

main ()

{
    int n_picks=125198, n_rec=5511, n_src=451;

    float *time, *offset; 
    int *source, *receiver,i,j,k,nelm;
    float *x_rec, *y_rec, *z_rec;
    float *x_src, *y_src, *z_src;
    float v0=1000, vref=5000, d0=300, d0s;
    float d0r, theta, dx, dy, dt, tij, offst;
    float *elem, t0;
    int *irow, *icol,rowindex;
    
    FILE *fp_pick, *fp_rec, *fp_src, *fp_dt, *fp_jac, *fp_src_rec_init;


/***get the picks
*/
    if( (fp_pick = fopen("picks_vconstant","r"))==NULL)
         fprintf(stderr,"can't open picks file 1 \n");

    time= (float *) calloc(sizeof(float),n_picks);
    offset= (float *) calloc(sizeof(float),n_picks);
    source= (int *) calloc(sizeof(float),n_picks);
    receiver= (int *) calloc(sizeof(float),n_picks);
			    

    fread( time, sizeof(float), n_picks, fp_pick);
    fread( offset,  sizeof(float), n_picks, fp_pick);
    fread( source, sizeof(int), n_picks, fp_pick);
    fread( receiver, sizeof(int), n_picks, fp_pick);

    fprintf(stderr,"Pick 1: time=%f, offset=%f, source=%d, receiver=%d\n",
                    time[0],offset[0],source[0],receiver[0]);



/***get the receiver (x,y,z)
*/

    if( (fp_rec = fopen("rec_xyz","r"))==NULL)
         fprintf(stderr,"can't open receiver file \n");

    x_rec= (float *) calloc(sizeof(float),n_rec);
    y_rec= (float *) calloc(sizeof(float),n_rec);
    z_rec= (float *) calloc(sizeof(float),n_rec);

    fread( x_rec, sizeof(float), n_rec, fp_rec);
    fread( y_rec, sizeof(float), n_rec, fp_rec);
    fread( z_rec, sizeof(float), n_rec, fp_rec);

    fprintf(stderr,"Receiver 1: x=%f, y=%f, z=%f\n",
                    x_rec[0],y_rec[0],z_rec[0]);



/***get the source (x,y,z)
*/

    if( (fp_src = fopen("src_xyz","r"))==NULL)
         fprintf(stderr,"can't open source file \n");

    x_src= (float *) calloc(sizeof(float),n_src);
    y_src= (float *) calloc(sizeof(float),n_src);
    z_src= (float *) calloc(sizeof(float),n_src);

    fread( x_src, sizeof(float), n_src, fp_src);
    fread( y_src, sizeof(float), n_src, fp_src);
    fread( z_src, sizeof(float), n_src, fp_src);

    fprintf(stderr,"Source 1: x=%f, y=%f, z=%f\n",
                    x_src[0],y_src[0],z_src[0]);


    if( (fp_dt = fopen("rhs.paul","w"))==NULL)
         fprintf(stderr,"can't open delta_t file 1 \n");
    if( (fp_src_rec_init = fopen("src_rec_init","w"))==NULL)
         fprintf(stderr,"can't open src_rec_init file 1 \n");

    nelm=0;
    theta = asin(v0/vref);
    for(i=0; i < n_picks; i++){
    if(offset[i] > 1000 && offset[i] < 2000){
/*	d0r = d0 - abs(z_rec[receiver[i]]);
	d0s = d0 - abs(z_src[source[i]]); */
        d0r = d0;
        d0s = d0;
	tij = 1/(v0 * cos(theta))*(d0s + d0r)
	    + (offset[i] - (d0s + d0r)*tan(theta))/vref;
	dt = time[i] - tij;
	nelm = nelm + 1;
	fprintf(fp_dt,"%f \n",dt);
    }
}

    for(i=0; i < n_src; i++){
	    fprintf(fp_src_rec_init,"%f \n",d0/v0);
	}
    for(i=0; i < n_rec; i++){
            fprintf(fp_src_rec_init,"%f \n",d0/v0);
        }

    nelm = 2 * nelm;
    fprintf(stdout,"number of elements = %d \n",nelm);
    elem = (float *) calloc(sizeof(float),nelm);
    irow = (int *) calloc(sizeof(float),nelm);
    icol = (int *) calloc(sizeof(float),nelm);

    k = 0;
    rowindex = 0;
    for(i=0; i < n_picks; i++){
	if(offset[i] > 1000 && offset[i] < 2000){
	elem[k] = 1;
	icol[k] = source[i]-1;
	irow[k] = rowindex;
	k=k+1;
	elem[k] = 1;
	icol[k] = receiver[i] + n_src -1;
	irow[k] = rowindex;
	k=k+1;
	rowindex = rowindex + 1;
/*	fprintf(stderr,"i = %d, k = %d \n",i,k);*/
    }
    }

    if( (fp_jac = fopen("matrix.paul","w"))==NULL)
         fprintf(stderr,"can't open matrix file \n");    

    for(i=0; i < nelm; i++) fprintf(fp_jac,"%f\n",elem[i]);
    for(i=0; i < nelm; i++) fprintf(fp_jac,"%d\n",irow[i]);
    for(i=0; i < nelm; i++) fprintf(fp_jac,"%d\n",icol[i]);
	



}