eval_phinc()

void eval_phinc	(	COMPLEX	cfreq,
		int	icounter
	)

This routine evaluates the total incident potential field at the control points and at microphones. If knw=2, it also evaluates the normal derivative of the total incident field at the boundary control points.

Parameters

[in]	cfreq	COMPLEX frequency
[in]	icounter	frequency index

                                           {

  int i;
  int npsiu,nphiu;
  int ipsiu,isourc,icntr,iplanw,imics;
  struct acoustic_point xs,ws;
  COMPLEX grgr;
  COMPLEX grg[3];
  DOUBLE rad,theta;
  DOUBLE wnum;
  DOUBLE espo;
  DOUBLE pi;
  COMPLEX zero,one;
  COMPLEX phincval,psincval,phincmval;

  int locRows;
  int gRows;

  int locMRows;
  int gMRows;
  int izero = 0;

  struct dipole dip;
  struct vector vdip;
  DOUBLE cosdelta;
  int idipole;

#ifdef _ACOUSTO_DEBUG
  FILE *f;
  char fname[MAX_PATH];
#endif

  nphiu = modgeominfo->ncntr;
  npsiu = modgeominfo->ncntr;

  gMRows = modgeominfo->nmics;
  gRows = nphiu;

  locRows  = numroc_(&gRows ,&runinfo->row_block_size,&runinfo->myrow,&izero,&runinfo->nprows);
  locMRows = numroc_(&gMRows,&runinfo->mics_block_size,&runinfo->myrow,&izero,&runinfo->nprows);

  if(icounter == 0){
   CREATE_VEC(solution->vecPhinc ,locRows);
   CREATE_VEC(solution->vecPsinc ,locRows);
   CREATE_VEC(solution->vecPhincm,locMRows);
  }

  VecZeroEntries(solution->vecPhinc ,locRows);
  VecZeroEntries(solution->vecPsinc ,locRows);
  VecZeroEntries(solution->vecPhincm,locMRows);

  if(runinfo->mycol != 0){ // vectors are distributed only on first col of grid
    return;
  }

  pi   = 4*atan(1.0);
  wnum = CIMAG(cfreq)/runinfo->vsound;
  zero = 0.0;
  one  = 1.0;



  for(i=0;i<locRows;i++){
    sc_l2g(i,runinfo->myrow,gRows,runinfo->nprows,runinfo->row_block_size,&icntr);
    ipsiu = icntr;

    /* sources contribution */
    psincval = 0.0 + I*0.0;
    phincval = 0.0 + I*0.0;
    for(isourc=0;isourc<modsolinfo->nsourc;isourc++){
      vec_copy(&xs.pos, solution->sources[isourc].pos);
      xs.amplitude = solution->sources[isourc].amplitude;
      delaysc2(geometry->cntr[icntr],xs.pos,&rad,&theta);

      phincval -= (xs.amplitude*CEXP(-cfreq*theta)/(4.*pi*rad));
      if(2 == modsolinfo->knw){
        grgr = -(xs.amplitude*(cfreq/runinfo->vsound+1.0/rad)*CEXP(-cfreq*theta)) / (4.*pi*rad*rad);

        grg[0] = grgr*(xs.pos.x - geometry->cntr[ipsiu].x);
        grg[1] = grgr*(xs.pos.y - geometry->cntr[ipsiu].y);
        grg[2] = grgr*(xs.pos.z - geometry->cntr[ipsiu].z);

        psincval += grg[0] * geometry->elements[ipsiu].n.x +
                    grg[1] * geometry->elements[ipsiu].n.y +
                    grg[2] * geometry->elements[ipsiu].n.z;

      }
    }
    /* dipoles contribution */
    for(idipole=0;idipole<modsolinfo->ndipoles;idipole++){
      vec_copy(&dip.pos, solution->dipoles[idipole].pos);
      vec_copy(&dip.dir, solution->dipoles[idipole].dir);
      dip.amplitude = solution->dipoles[idipole].amplitude;
      dip.d         = solution->dipoles[idipole].d;
      delaysc2(geometry->cntr[icntr],dip.pos,&rad,&theta);

      vec_diff(&vdip, geometry->cntr[icntr],dip.pos);
      cosdelta = vec_dot(dip.dir,vdip)/(vec_mod(vdip));
      LOG_DIPOLE(LOG_DEBUG,idipole,dip);
      phincval += (dip.amplitude*dip.d)*cosdelta*((cfreq/(runinfo->vsound*rad)) - (1.0/(rad*rad)))* CEXP(-cfreq*theta) /(4*pi);
//      phincval += (dip.amplitude*dip.d)*cosdelta*((cfreq/(runinfo->vsound*rad)))* CEXP(-cfreq*theta) /(4*pi);

      if(2 == modsolinfo->knw){
        grgr = ((dip.amplitude*dip.d*cosdelta)/(4*pi))*((2.0/(rad*rad*rad)) - ((cfreq*cfreq)/(runinfo->vsound*runinfo->vsound)))*CEXP(-cfreq*theta);

        grg[0] = grgr*(xs.pos.x - geometry->cntr[ipsiu].x);
        grg[1] = grgr*(xs.pos.y - geometry->cntr[ipsiu].y);
        grg[2] = grgr*(xs.pos.z - geometry->cntr[ipsiu].z);

        psincval += grg[0] * geometry->elements[ipsiu].n.x +
                    grg[1] * geometry->elements[ipsiu].n.y +
                    grg[2] * geometry->elements[ipsiu].n.z;

      }

    }

    /* planar waves contribution */
    for(iplanw=0;iplanw<modsolinfo->nplanw;iplanw++){
      vec_copy(&ws.pos,solution->planw[iplanw].pos);
      vec_scale(&ws.pos,wnum);
      ws.amplitude = solution->planw[iplanw].amplitude;
      espo = vec_dot(ws.pos,geometry->cntr[icntr]);

      phincval += ws.amplitude*CEXP(-I*espo);


      if(2 == modsolinfo->knw){
        grgr = - I * ws.amplitude * CEXP(-I*espo);

        grg[0] = grgr * ws.pos.x;
        grg[1] = grgr * ws.pos.y;
        grg[2] = grgr * ws.pos.z;

        psincval +=  grg[0] * geometry->elements[ipsiu].n.x +
                     grg[1] * geometry->elements[ipsiu].n.y +
                     grg[2] * geometry->elements[ipsiu].n.z;
      }

    }
    solution->vecPsinc[i] = psincval;
    solution->vecPhinc[i] = phincval;

  }

  /* microphones */
  for(i=0;i<locMRows;i++){
    sc_l2g(i,runinfo->myrow,gMRows,runinfo->nprows,runinfo->mics_block_size,&imics);
    phincmval = 0.0 +I*0.0;
    // sources
    for(isourc=0;isourc<modsolinfo->nsourc;isourc++){
      vec_copy(&xs.pos,solution->sources[isourc].pos);
      xs.amplitude = solution->sources[isourc].amplitude;
      delaysc2(geometry->mics[imics],xs.pos,&rad,&theta);
      phincmval -= xs.amplitude*CEXP(-cfreq*theta) / (4.0*pi*rad);
    }
    // dipoles
    for(idipole=0;idipole<modsolinfo->ndipoles;idipole++){
      vec_copy(&dip.pos, solution->dipoles[idipole].pos);
      vec_copy(&dip.dir, solution->dipoles[idipole].dir);
      dip.amplitude = solution->dipoles[idipole].amplitude;
      dip.d         = solution->dipoles[idipole].d;
      delaysc2(geometry->mics[imics],dip.pos,&rad,&theta);

      vec_diff(&vdip, geometry->mics[imics],dip.pos);
      cosdelta = vec_dot(dip.dir,vdip)/(vec_mod(vdip));
      phincmval += (dip.amplitude * dip.d)*cosdelta*((cfreq/(runinfo->vsound*rad)) - (1.0/(rad*rad)))* CEXP(-cfreq*theta) /(4*pi);
//      phincmval += (dip.amplitude*dip.d)*cosdelta*((cfreq/(runinfo->vsound*rad)))* CEXP(-cfreq*theta) /(4*pi);

    }
    // planar waves
    for(iplanw=0;iplanw<modsolinfo->nplanw;iplanw++){
      vec_copy(&ws.pos,solution->planw[iplanw].pos);
      vec_scale(&ws.pos,wnum);
      ws.amplitude = solution->planw[iplanw].amplitude;
      espo = vec_dot(ws.pos,geometry->mics[imics]);
      phincmval += ws.amplitude*CEXP(-I*espo);
    }
    solution->vecPhincm[i] = phincmval;
  }


//#define _ACOUSTO_DEBUG
#ifdef _ACOUSTO_DEBUG
  FILE *f;
  char fname[MAX_PATH];
  sprintf(fname,"psi.%d.out",modgeominfo->nelmb);
  f = fopen(fname,"w");
  for(i=0;i<locRows;i++){
    sc_l2g(i,runinfo->myrow,gRows,runinfo->nprows,runinfo->row_block_size,&icntr);
    fprintf(f," %d %f %f %f %f %f %f %f \n",icntr,
                           geometry->cntr[icntr].x,
                           geometry->cntr[icntr].y,
                           geometry->cntr[icntr].z,
                           CREAL(solution->vecPhinc[i]),CIMAG(solution->vecPhinc[i]),
                           CREAL(solution->vecPsinc[i]),CIMAG(solution->vecPsinc[i])
                           );
  }
  fclose(f);
  

  sprintf(fname,"phincm.%d%d",runinfo->myrow,runinfo->mycol);
  f = fopen(fname,"w");
  for(i=0;i<locMRows;i++){
    sc_l2g(i,runinfo->myrow,gMRows,runinfo->nprows,runinfo->mics_block_size,&imics);
    fprintf(f," %d %f %f  \n",imics,
                           CREAL(solution->vecPhincm[i]),CIMAG(solution->vecPhincm[i])
                           );
  }
  fclose(f);
#endif
//#undef _ACOUSTO_DEBUG
  
 
}