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

                         SYRTHES version 5.0
                         -------------------

     This file is part of the SYRTHES Kernel, element of the
     thermal code SYRTHES.

     Copyright (C) 2009 EDF S.A., France

     contact: syrthes-support@edf.fr


     The SYRTHES Kernel 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 3 of
     the License, or (at your option) any later version.

     The SYRTHES Kernel 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.


     You should have received a copy of the GNU General Public License
     along with the SYRTHES Kernel; if not, write to the
     Free Software Foundation, Inc.,
     51 Franklin St, Fifth Floor,
     Boston, MA  02110-1301  USA

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

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

#include "syr_usertype.h"
# include "syr_tree.h"
# include "syr_abs.h"
# include "syr_meteo.h"
# include "syr_option.h"
# include "syr_const.h"
# include "syr_bd.h"
# include "syr_hmt_bd.h"
# include "syr_parall.h"
#include "syr_fluid0d.h"
#include "syr_fluid1d.h"
#include "syr_couplageSR.h"
#include "syr_soleil.h"
#include "syr_bilan.h"

# include "syr_resray.h"
#include "syr_utilitaire.h"
#include "syr_operations.h"
#include "syr_generic.h"
#include "syr_user_fct.h"
#include "syr_user_ray.h"

#ifdef _SYRTHES_MPI_
# include <mpi.h>
static MPI_Status status;
#endif

extern struct Affichages affich;
extern struct Performances perfo;
extern int lray_fdf;
extern int lfluid0d;

double epsgcs_rad=1.e-6;
int optim=0;

int ndecoup_max=0;
extern FILE  *fchror;
extern char nomresucr[CHLONG];

/*|======================================================================|
  | SYRTHES 5.0                                       COPYRIGHT EDF 2009 |
  |======================================================================|
  | AUTEURS  : C. PENIGUEL, I. RUPP                                      |
  |======================================================================|
  | radiation                                                            |
  |                                                                      |
  |======================================================================| */
void radiation(struct PasDeTemps *pasdetemps,struct GestionFichiers gestionfichiers,
               struct Maillage maillnodes,struct Maillage maillnodray,
               struct MaillageCL maillnodeus,
               double *tmps,double *tmpsa,struct Couple scoupr,struct Couple rcoups,
               struct FacForme ***listefdf,double *diagfdf,struct FacGebhart ***listefdg,double **diagfdg,
               double *tmpray,double **radios,double **firay,double *trayeq,double *erayeq,
               struct ProphyRay phyray,
               double **epropr,struct Clim fimpray,
               struct PropInfini *propinf,struct MatriceRay matrray,
               struct Soleil *soleil,struct Bilan bilansolaire,
               struct Meteo meteo,struct Myfile myfile,
               struct Compconnexe volconnexe,struct Vitre vitre,
               struct Fluid0d fluid0d,struct SymPer raysymper,
               struct Travail trav,
               struct SDparall sdparall,struct SDparall_ray sdparall_ray)
{
  int i;
  double tcpu1,tcpu2;
  double *trav1;
  char chi[5],chn[20];


  pasdetemps->nbparay++;

  /* isa ?????????????????????????? */
/*   trav1=trav.tab[0]; */
/*   for (i=0;i<maillnodes.npoin;i++)  *(trav1+i)=*(tmps[ADR_T]+i) */
/*     +0.5*pasdetemps->rdtts/pasdetemps->rdttsprec*(*(tmps[ADR_T]+i)-*(tmpsa[ADR_T]+i)); */
/*   if (sdparall_ray.npartsray>1) */
/*     transSR_parall(maillnodes.ndim,maillnodeus.node,rcoups,trav1,tmpray,sdparall_ray); */
/*   else */
/*     transSR(maillnodes.ndim,maillnodeus.node,rcoups,trav1,tmpray); */


  /* le transfert se fait avec tous les proc */
  /* --------------------------------------- */
  if (sdparall_ray.nparts>1)
    transSR_parall(maillnodes.ndim,maillnodeus.node,rcoups,tmps,tmpray,sdparall_ray);
  else
    transSR(maillnodes.ndim,maillnodeus.node,rcoups,tmps,tmpray);



  /* la resolution ne se fait que sur les noeuds de ray */
  /* -------------------------------------------------- */
  if (sdparall_ray.rangray>-1)
    {

      user_ray(maillnodray,tmpray,phyray,propinf,fimpray,vitre,*pasdetemps);
      verif_ray(maillnodray,phyray);


      /* mise à jour des flux solaires */
       if (soleil->actif) resol_soleil(soleil,maillnodray,phyray,
                                       *pasdetemps,vitre,meteo,myfile,*propinf,sdparall);


      /* resolution du rayonnement */
      tcpu1=cpusyrt();

      if (lray_fdf==1)
        resray(pasdetemps->ntsyr,maillnodray,listefdf,diagfdf,tmpray,
               radios,firay,trayeq,erayeq,phyray,epropr,
               fimpray,*propinf,
               matrray,*soleil,bilansolaire,
               volconnexe,vitre,trav.tab[0],sdparall_ray);
      else
        {
          for (i=0;i<maillnodray.nelem;i++) *(tmpray+i) += tkel;

          if (sdparall_ray.npartsray==1)
            fi2teqfdg(maillnodray.nelem,tmpray,listefdg[0],diagfdg,
                      maillnodray.volume,phyray,firay,trayeq,erayeq,
                      *propinf,fimpray,*soleil,vitre,fluid0d,raysymper);
          else
            fi2teqfdg_parall(maillnodray.nelem,tmpray,listefdg,diagfdg,
                             maillnodray.volume,phyray,firay,trayeq,erayeq,
                             *propinf,fimpray,*soleil,vitre,fluid0d,raysymper,
                             trav,sdparall_ray);
      printf("ggggggggggg T=%f\n",tmpray[33]);

          for (i=0;i<maillnodray.nelem;i++)
            {*(tmpray+i) -= tkel; *(trayeq+i) -= tkel;}

        }

      /* bilans sur les flux solaires */
      if (soleil->actif && bilansolaire.nb)
        bilsolaire(maillnodray,*soleil,bilansolaire,phyray.bandespec,propinf->fdfnp1,pasdetemps->tempss);

      tcpu2=cpusyrt();
      perfo.cpu_1iter_ray=tcpu2-tcpu1;
      perfo.cpu_iter_ray+=perfo.cpu_1iter_ray;


    }

  /* le transfert se fait avec tous les proc */
  /* --------------------------------------- */
  if (sdparall_ray.nparts>1)
    transRS_parall(maillnodeus.ndmat,scoupr,trayeq,erayeq,sdparall_ray);
  else
    transRS(maillnodeus.ndmat,scoupr,trayeq,erayeq);


}

/*|======================================================================|
  | SYRTHES 5.0                                       COPYRIGHT EDF 2009 |
  |======================================================================|
  | AUTEURS  : C. PENIGUEL, I. RUPP                                      |
  |======================================================================|
  | resray                                                               |
  |                                                                      |
  |======================================================================| */
void resray(int ntsyr,struct Maillage maillnodray,
            struct FacForme ***listefdf,double *diagfdf,double *tmpray,
            double **radios,double **firay,double *trayeq,double *erayeq,
            struct ProphyRay phyray,
            double **epropr,struct Clim fimpray,
            struct PropInfini propinf,struct MatriceRay matrray,
            struct Soleil soleil,struct Bilan bilansolaire,
            struct Compconnexe volconnexe,struct Vitre vitre,
            double *trav0,struct SDparall_ray sdparall_ray)
{

  int numbs,n,ngfac,i,j,niter,nel,nume,iv,nf;
  double v,*r;
  static int lprem=1;
  double c2=1.4388e-2,x1,x2,w1,w2,xx;

  for (i=0;i<maillnodray.nelem;i++) *(tmpray+i) += tkel;

  if (lprem)
    {
      lprem=0;
      if (phyray.bandespec.nb==1)
        for (i=0;i<maillnodray.nelem;i++)
          radios[0][i]=sigma*phyray.emissi[0][i]*pow(tmpray[i],4);
      else
        for (numbs=0;numbs<phyray.bandespec.nb;numbs++)
          {
            x1=c2/phyray.bandespec.borneinf[numbs]; x2=c2/phyray.bandespec.bornesup[numbs];
            for (i=0;i<maillnodray.nelem;i++)
              {
                v=x1/tmpray[i]; w1=wiebel(v);
                v=x2/tmpray[i]; w2=wiebel(v);
                radios[numbs][i]=sigma*phyray.emissi[numbs][i]*pow(tmpray[i],4)*(w2-w1);
              }
          }
    }

  /* second membre */
  /* ------------- */
  for (numbs=0;numbs<phyray.bandespec.nb;numbs++)
    {
      smbray(numbs,maillnodray.nelem,maillnodray.volume,tmpray,propinf,phyray,
             epropr,radios,vitre);

      for (i=0;i<fimpray.nelem;i++)
        epropr[numbs][fimpray.numf[i]]=fimpray.val1[numbs][i]*maillnodray.volume[fimpray.numf[i]];


      if (soleil.actif && (soleil.model==3 || soleil.h*180/Pi>=1.))
        {
          /* flux solaire diffus sur toutes les facettes de la voute */
          for (i=0;i<soleil.voute->nelem;i++)
            {
              nf=soleil.voute->numf[i];
              epropr[numbs][nf] = soleil.voute->fluxdiffus[numbs] * maillnodray.volume[nf];
            }

          /* flux solaire direct surajouté sur LA facette correspondante au  soleil */
          nf=soleil.voute->numf[soleil.voute->numFacePosSoleil];
          epropr[numbs][nf] += soleil.voute->fluxdirect[numbs] * maillnodray.volume[nf];

        }
    }

  /* resolution du systeme */
  /* --------------------- */
  if (vitre.actif)
    {
      if (SYRTHES_LANG == FR)
        printf("Traitement des vitres non disponible\n");
      else if (SYRTHES_LANG == EN)
        printf("Windows handling not yet completly implemented\n");
      syrthes_exit(1);

    }
  else
    {
      for (j=0;j<phyray.bandespec.nb;j++)
        {
          if (syrglob_nparts==1 || syrglob_rang==0) {
            if (SYRTHES_LANG == FR) printf("       Bande spectrale %d\n",j+1);
            else if (SYRTHES_LANG == EN) printf("       Spectral band %d\n",j+1);
          }
          for (i=0;i<fimpray.nelem;i++)
            {phyray.reflec[j][fimpray.numf[i]]=1.;
              phyray.absorp[j][fimpray.numf[i]]=0.;}


          if (soleil.actif)
            {
              for (i=0;i<soleil.voute->nelem;i++)
                {phyray.reflec[j][soleil.voute->numf[i]]=1.;
                  phyray.absorp[j][soleil.voute->numf[i]]=0.;}
            }

          if (sdparall_ray.npartsray>1)
            rrayrc_parall(maillnodray.nelem,listefdf,diagfdf,
                          maillnodray.volume,phyray.reflec[j],epropr[j],
                          matrray,&niter,trav0,sdparall_ray);
          else
            rrayrc(maillnodray.nelem,listefdf[0],diagfdf,
                   maillnodray.volume,phyray.reflec[j],epropr[j],
                   matrray,&niter);

          if (niter) for (i=0;i<maillnodray.nelem;i++) radios[j][i]=matrray.x[i];
        }
    }


  if (affich.ray_resray)
    {
      if (SYRTHES_LANG == FR)
        printf("Fin de rrayrc : impression de la radiosite pour la bande 0\n");
      else if (SYRTHES_LANG == EN)
        printf("End of rrayrc : printing for radiosity of band 0\n");
      for (i=0;i<maillnodray.nelem;i++) printf("radios[%d]=%25.18e\n",i,radios[0][i]);
    }

  /* 3- preparation des donnees equivalentes */
  /* --------------------------------------- */
  if (sdparall_ray.npartsray>1)
    /* b et di et xm1 sont des vect auxiliaires pour rrayrc : */
    /* on s'en sert ici comme tableaux de travail ! */
    fi2teq_parall(maillnodray.nelem,tmpray,listefdf,diagfdf,
                  maillnodray.volume,phyray,radios,
                  firay,trayeq,erayeq,propinf,fimpray,soleil,vitre,
                  sdparall_ray,trav0,matrray.di,matrray.xm1);
  else
    /* b et di sont des vect auxiliaires pour rrayrc : on s'en sert ici comme tableaux de travail */
    fi2teq(maillnodray.nelem,tmpray,listefdf[0],diagfdf,
           maillnodray.volume,phyray,radios,
           firay,trayeq,erayeq,propinf,fimpray,soleil,vitre);


  for (i=0;i<maillnodray.nelem;i++)
    {*(tmpray+i) -= tkel; *(trayeq+i) -= tkel;}

}

/*|======================================================================|
  | SYRTHES 5.0                                       COPYRIGHT EDF 2009 |
  |======================================================================|
  | AUTEURS  : C. PENIGUEL, I. RUPP                                      |
  |======================================================================|
  | smbray                                                               |
  |                                                                      |
  |======================================================================| */
void smbray(int numbs,int nelray,
            double *sufray,double *tmpray,struct PropInfini propinf,
            struct ProphyRay phyray,double **epropr,double **radios,
            struct Vitre vitre)
{
  int iv,n,i,ng,ngvois,j,nel,nume;
  double c2,x1,x2,v,w1,w2,xj,*e,*r,x;
  double *fdfnp1;

  fdfnp1=propinf.fdfnp1;

  c2=1.4388e-2;
  if (propinf.actif||vitre.actif)
    {
      e=(double*)malloc(nelray*sizeof(double)); verif_alloue_double1d("smbray",e);
      r=(double*)malloc(nelray*sizeof(double)); verif_alloue_double1d("smbray",r);
      for (i=0;i<nelray;i++) {e[i]=0;r[i]=0;}
    }


  /* 1- cas d'une seule bande spectrale*/
  if (phyray.bandespec.nb==1)
    {
      for (i=0;i<nelray;i++)
        epropr[0][i]=sigma*phyray.emissi[0][i]*pow(tmpray[i],4)*sufray[i];


      if (propinf.actif)
        for (i=0;i<nelray;i++) e[i]+=fdfnp1[i]*sigma * propinf.EO[0][i] * pow(propinf.TO[i]+tkel,4);

    }

  /* 2- cas de bandes spectrales multiples*/
  else
    {
      x1=c2/phyray.bandespec.borneinf[numbs]; x2=c2/phyray.bandespec.bornesup[numbs];
      for (i=0;i<nelray;i++)
        {
          v=x1/tmpray[i]; w1=wiebel(v);
          v=x2/tmpray[i]; w2=wiebel(v);
          epropr[numbs][i]=sigma*phyray.emissi[numbs][i]*pow(tmpray[i],4)*(w2-w1)*sufray[i];
        }
      if (propinf.actif)
        for (i=0;i<nelray;i++)
          {
            v=x1/(propinf.TO[i]+tkel);  w1=wiebel(v);
            v=x2/(propinf.TO[i]+tkel);  w2=wiebel(v);
            e[i]+=fdfnp1[i]*sigma * propinf.EO[numbs][i] * pow(propinf.TO[i]+tkel,4)*(w2-w1);
          }

    }

  if (propinf.actif)
    for (i=0;i<nelray;i++)
      epropr[numbs][i]+=phyray.reflec[numbs][i]*e[i];

  if (vitre.actif)
    {
      for (i=0;i<nelray;i++)
        if ((iv=vitre.voisin[i])>-1)
          epropr[numbs][i]+=phyray.transm[numbs][i]*e[iv];
    }

  /*   printf("SMBRAY second membre epropr\n");
       for (i=0;i<nelray;i++) printf(" i=%d epropr/S=%f\n",i,epropr[0][i]/sufray[i]); */

  if (propinf.actif||vitre.actif) {free(e); free(r);}

}
/*|======================================================================|
  | SYRTHES 5.0                                       COPYRIGHT EDF 2009 |
  |======================================================================|
  | AUTEURS  : C. PENIGUEL, I. RUPP                                      |
  |======================================================================|
  | wiebel                                                               |
  |                                                                      |
  |======================================================================| */
double wiebel(double v)
{
  int m;
  double w,v2,v4;
  const double pi4=.153989717364e+00;
  const double z0=.333333343267e+00;
  const double z1=.125000000000e+00;
  const double z2=.166666675359e-01;
  const double z4=.198412701138e-03;
  const double z6=.367430925508e-05;
  const double z8=.751563220547e-07;

  if (v>=2.)
    {
      w=0;
      for (m=1;m<6;m++) w+=exp(-m*v)/(m*m*m*m) * (((m*v+3)*m*v+6)*m*v+6);
      w=w * pi4;
    }
  else
    {
      v2=v*v;
      v4=v2*v2;
      w=z0 - z1*v + z2*v2 - z4*v2*v2 + z6*v4*v2 - z8*v4*v4;
      w=1. - pi4*v2*v*w;
    }
  return w;
}

/*|======================================================================|
  | SYRTHES 5.0                                       COPYRIGHT EDF 2009 |
  |======================================================================|
  | AUTEURS  : C. PENIGUEL, I. RUPP                                      |
  |======================================================================|
  | gausei                                                               |
  |                                                                      |
  |======================================================================| */
void rrayrc(int nelray,struct FacForme **listefdf,double *diagfdf,
            double *sufray,
            double *rho,double *epr,
            struct MatriceRay matrray,int *niter)
{
  int n,ngfac,nn;
  double aa,bb,cc,d,ee,alfa,aalfa,tau,alfam1,alftau;
  double aux1,aux2,aux,denom;
  int i,j,nitsmo;
  double x0,resnor,prsca1;
  double epsis,zero,epssmo;
  double *xx;
  struct FacForme *pff;
  double *x,*b,*xm1,*gd,*res,*z,*di,*resm1;


  x     = matrray.x;
  b     = matrray.b;
  xm1   = matrray.xm1;
  gd    = matrray.gd;
  res   = matrray.res;
  z     = matrray.z;
  di    = matrray.di;
  resm1 = matrray.resm1;

  zero=0.;
  nitsmo=100;
  epssmo=1.e-12;
  n=0;

  /* 1- initialisation des vecteurs auxiliaires */
  for (i=0;i<nelray;i++)
    {
      xm1[i]=x[i]=epr[i];
      b[i]=epr[i];
      di[i]=sufray[i]-rho[i]*diagfdf[i];
    }

  /* norme du second membre */
  for (i=0,prsca1=0,xx=x;i<nelray;i++,xx++) prsca1+= *xx * *xx;
  x0=sqrt(prsca1);

  if (x0< 1.e-20 ) x0=1.e-6;
  epsis=1.e-4 * x0;

  for (i=0;i<nelray;i++)
    {
      res[i]=0.;
      for (nn=0;nn<listefdf[i]->nb;nn++)
        res[i]-=listefdf[i]->fdf[nn]*rho[i]*x[listefdf[i]->numenface[nn]];

      res[i] += di[i]*x[i]-b[i];
   }


  for (i=0,prsca1=0,xx=res;i<nelray;i++,xx++) prsca1+= *xx * *xx;
  resnor=sqrt(prsca1);
/*   printf(" >>>> ray=resnor_ini seq : %25.18e \n",resnor); */

  if (SYRTHES_LANG == FR)
    printf("       RRAYRC : Iteration   Precision relative  Precision absolue\n");
  else if (SYRTHES_LANG == EN)
    printf("       RRAYRC : Iteration   Relative Precision  Absolute Precision\n");


  if ( resnor<=epsis && resnor<=epsgcs_rad*sqrt((double)(nelray)))
    {
      printf("                 %4d         %12.5e                %12.5e\n",
             n,resnor/x0,resnor/sqrt((double)(nelray)));
      *niter=n;
    }

  else
    {
      for (i=0;i<nelray;i++) resm1[i]=res[i];
      n=0;

      do
        {
          n++;
          for (i=0;i<nelray;i++) gd[i]=res[i]/di[i];
          for (i=0;i<nelray;i++)
            {
              z[i]=0.;
              for (nn=0;nn<listefdf[i]->nb;nn++) z[i]-=listefdf[i]->fdf[nn]*rho[i]*gd[listefdf[i]->numenface[nn]];
              z[i] += di[i]*gd[i];
            }
          aa=0.;bb=0.;cc=0.;d =0.;ee=0.;
          for (i=0;i<nelray;i++)
            {
              aa += z[i]*res[i];
              bb += z[i]*resm1[i];
              aux=res[i]-resm1[i];
              cc += res[i]*aux;
              d  += resm1[i]*aux;
              ee +=z[i]*z[i];
            }

          denom=(cc-d )*ee-(aa-bb)*(aa-bb);
          if(fabs(denom)< 1.e-20) alfa=1.;
          else alfa=((aa-bb)*bb-d *ee)/denom;
          aalfa=fabs(alfa);
          if (aalfa<1e-20 || fabs(aalfa-1.)< 1.e-20) {alfa=1.; tau=aa/ee;}
          else tau=aa/ee + (1.-alfa)/alfa * bb/ee;
          alfam1 =1.-alfa;
          alftau =-alfa*tau;
          for (i=0;i<nelray;i++)
            {
              aux1    =res[i];
              aux2    =x[i];
              res[i]  =alfa*aux1+alfam1*resm1[i]+alftau*z[i];
              resm1[i]=aux1;
              x[i]    =alfa*aux2+alfam1*xm1[i]+alftau*gd[i];
              xm1[i]  =aux2;
            }

          for (i=0,prsca1=0,xx=res;i<nelray;i++,xx++) prsca1+= *xx * *xx;
          resnor=sqrt(prsca1);
/*        printf("rrayrc_seq : it=%d   resnor=%25.18e \n",n,resnor); */

          if (n%10==0)
            printf("                 %4d         %12.5e        %12.5e\n",
                   n,resnor/x0,resnor/sqrt((double)(nelray)));

        }
      while(!( (resnor<=epsis &&  resnor<=epsgcs_rad*sqrt((double)(nelray))) || n>=nitsmo));

      if (n%10!=0)
        printf("                 %4d         %12.5e        %12.5e\n",
               n,resnor/x0,resnor/sqrt((double)(nelray)));
      *niter=n;
    }
/*   for (i=0;i<nelray;i++)  */
/*     printf("fin de rrayrc proc=%d temray[%d]=%f\n",syrglob_rang,i,matrray.x[i]); */
}



/*|======================================================================|
  | SYRTHES 5.0                                       COPYRIGHT EDF 2009 |
  |======================================================================|
  | AUTEURS  : C. PENIGUEL, I. RUPP                                      |
  |======================================================================|
  | fi2teq                                                               |
  |                                                                      |
  |======================================================================| */
void fi2teq(int nelray,double *tmpray,
            struct FacForme **listefdf,double *diagfdf,
            double *sufray,
            struct ProphyRay phyray,
            double **radios,double **firay,
            double *trayeq,double *erayeq,
            struct PropInfini propinf,struct Clim fimpray,
            struct Soleil soleil,
            struct Vitre vitre)
{
  int nb,nn,n,m,i,j,ngfac,nel,nume,iv,jv,numbs,nf;
  double eps;
  double w1,w2,tfac,xj1;
  double ee,c2,x1,x2,x,v,v2,v4,rr,**e,**ev,xx,f;
  struct FacForme *pff;
  double *fdfnp1;

  fdfnp1=propinf.fdfnp1;

  if (phyray.bandespec.nb>100)
    {
      if (SYRTHES_LANG == FR)
        {
          printf("fi2teq : le nombre de bandes spectrales ne peut etre superieur a 100\n");
          printf("         dans le cas improbable ou cela pose pb contacter CP/IR \n");
        }
      else if (SYRTHES_LANG == EN)
        {
          printf("fi2teq : the number of spectral band cannot be above 100\n");
          printf("         in the very unlikely case it induces a pb contact  CP/IR \n");
        }
      syrthes_exit(1);
    }

  eps=1e-6;
  c2=1.4388e-2;
  for (i=0;i<nelray;i++) trayeq[i]=erayeq[i]=0.;

  e=(double**)malloc(phyray.bandespec.nb*sizeof(double*));
  for (n=0;n<phyray.bandespec.nb;n++) e[n]=(double*)malloc(nelray*sizeof(double));
  verif_alloue_double2d(phyray.bandespec.nb,"fi2teq",e);
  for (i=0;i<nelray;i++) for (n=0;n<phyray.bandespec.nb;n++) {e[n][i]=0;}


  /* emissivite equivalente */
  if (phyray.bandespec.nb==1)
    {
      for (i=0;i<nelray;i++) erayeq[i]=phyray.emissi[0][i];
    }
  else
    for (n=0;n<phyray.bandespec.nb;n++)
      {
        x1=c2/phyray.bandespec.borneinf[n]; x2=c2/phyray.bandespec.bornesup[n];

        for (i=0;i<nelray;i++)
          {
            v=x1/tmpray[i]; w1=wiebel(v);
            v=x2/tmpray[i]; w2=wiebel(v);
            erayeq[i] += (w2-w1)*phyray.emissi[n][i];
          }
      }


  /* temperature equivalente */
  if (phyray.bandespec.nb==1)
    {
      for (i=0;i<nelray;i++)
        {
          rr=0;
          for (nn=0;nn<listefdf[i]->nb;nn++)
            rr+=listefdf[i]->fdf[nn]*radios[0][listefdf[i]->numenface[nn]];

          rr+=diagfdf[i]*radios[0][i];

          if (propinf.actif)
            rr += fdfnp1[i]*sigma*propinf.EO[0][i] * pow(propinf.TO[i]+tkel,4);

          if (rr < 0.0) rr = 0.0;
          trayeq[i]=pow(rr/(sufray[i]*sigma),0.25);
        }
    }
  else
    {
      for (i=0;i<nelray;i++)
        {
          rr=0.;
          for (nb=0;nb<phyray.bandespec.nb;nb++)
            {
              xx=0;
              for (nn=0;nn<listefdf[i]->nb;nn++)
                xx+=listefdf[i]->fdf[nn]*radios[nb][listefdf[i]->numenface[nn]];

              xx+=diagfdf[i]*radios[nb][i];
              if (propinf.actif)
                {
                  x1=c2/phyray.bandespec.borneinf[nb]; x2=c2/phyray.bandespec.bornesup[nb];

                  v=x1/(propinf.TO[i]+tkel);  w1=wiebel(v);
                  v=x2/(propinf.TO[i]+tkel);  w2=wiebel(v);

                  xx += fdfnp1[i]*sigma*propinf.EO[nb][i]*pow(propinf.TO[i]+tkel,4)*(w2-w1);
                }
              rr+=xx*phyray.emissi[nb][i];
            }
          trayeq[i]=pow(rr/(sufray[i]*erayeq[i]*sigma),0.25);
        }
    }



  /* 4.- calcul du flux de rayonnemnent (pour le post-processing) */
  for (nb=0;nb<phyray.bandespec.nb;nb++)
    {
      for (i=0;i<nelray;i++)
        {
          rr=0;
          for (nn=0;nn<listefdf[i]->nb;nn++)
            rr+=listefdf[i]->fdf[nn]*radios[nb][listefdf[i]->numenface[nn]];

          rr+=diagfdf[i]*radios[nb][i];

          if (propinf.actif)
            {
              x1=c2/phyray.bandespec.borneinf[nb]; x2=c2/phyray.bandespec.bornesup[nb];
              v=x1/(propinf.TO[i]+tkel);  w1=wiebel(v);
              v=x2/(propinf.TO[i]+tkel);  w2=wiebel(v);

              rr += fdfnp1[i]*sigma*propinf.EO[nb][i]*pow(propinf.TO[i]+tkel,4)*(w2-w1);
            }

          if (phyray.bandespec.nb==1)
            firay[nb][i]= phyray.emissi[nb][i]*(sigma*pow(tmpray[i],4)-rr/sufray[i]);
          else
            {
              v=x1/tmpray[i]; w1=wiebel(v);
              v=x2/tmpray[i]; w2=wiebel(v);
              firay[nb][i]= phyray.emissi[nb][i]*(sigma*(w2-w1)*pow(tmpray[i],4)-rr/sufray[i]);
            }
        }
    }

  /*
    if (vitre.actif)
    if (bandespec.nb==1)
    firay[0][i]= emissi[0][i]*sigma*pow(tmpray[i],4) - absorp[0][i]*ev[0][i];
    else
    for (nb=0;nb<bandespec.nb;nb++)
    for (i=0;i<nelray;i++)
    {
    v=x1/tmpray[ig]; w1=wiebel(v);
    v=x2/tmpray[ig]; w2=wiebel(v);
    firay[nb][i]= emissi[nb][i]*sigma*pow(tmpray[i],4) - absorp[n][i]*ev[n][i];
    }
    else
    if (bandespec.nb==1)
    firay[0][i]= emissi[0][i]*sigma*pow(tmpray[i],4) - absorp[0][i]*e[0][i];
    else
    for (nb=0;nb<bandespec.nb;nb++)
    for (i=0;i<nelray;i++)
    {
    v=x1/tmpray[ig]; w1=wiebel(v);
    v=x2/tmpray[ig]; w2=wiebel(v);
    firay[nb][i]= emissi[nb][i]*sigma*pow(tmpray[i],4) - absorp[nb][i]*e[nb][i];
    }
    */



  /* 5.1-mise a jour des facettes avec flux imposee */
  for (n=0;n<phyray.bandespec.nb;n++)
    for (i=0;i<fimpray.nelem;i++)
      {
        ngfac=fimpray.numf[i];
        firay[n][ngfac]=fimpray.val1[n][i];
      }

  /* 5.2-calcul de la temperature de la facette necessaire */
  for (i=0;i<fimpray.nelem;i++)
    {
      tfac=0.;
      ngfac=fimpray.numf[i];
      for (n=0;n<phyray.bandespec.nb;n++)
        tfac += (1.-phyray.emissi[n][ngfac])/phyray.emissi[n][ngfac]*fimpray.val1[n][i] + radios[n][ngfac]; // ATTENTION isa : il y avait i ici
      tmpray[ngfac]=pow(tfac/sigma,0.25);
      //printf("flux imp i=%d tmpray[%d]=%f radios=%f\n",i, ngfac,tmpray[ngfac], radios[0][ngfac]);
    }


  /* 5.2-calcul de la temperature de la facette necessaire pour la voute */
  double xflux;
  if (soleil.actif && (soleil.model==3 || soleil.h*180/Pi>=1.))
    {
      for (i=0;i<soleil.voute->nelem;i++)
        {
          tfac=0.;
          nf=soleil.voute->numf[i];
          for (n=0;n<phyray.bandespec.nb;n++)
            {
              xflux=soleil.voute->fluxdiffus[n];
              if (nf== soleil.voute->numf[soleil.voute->numFacePosSoleil]) xflux+=soleil.voute->fluxdirect[n];

              tfac += (1.-phyray.emissi[n][nf])/phyray.emissi[n][nf]*xflux +  radios[n][nf];
            }
          tmpray[nf]=pow(tfac/sigma,0.25);
          //printf("voute i=%d tmpray[%d]=%f radios=%f\n",i, nf,tmpray[nf],radios[0][ngfac]);

        }
    }




  if (affich.ray_fi2teq)
    {
      if (SYRTHES_LANG == FR)
        {
          printf("\n *** fi2teq : calcul des emissivites equivalentes ");
          printf(" des temperatures de rayonnement equivalentes\n");
          printf("             facette    emissivite equi   temp equi (degres c)\n");
        }
      else if (SYRTHES_LANG == EN)
        {
          printf("\n *** fi2teq : calculation of equivalent emissivity ");
          printf(" and of equivalent radiation temperature\n");
          printf("             face    emissivity equi   temp equi (degres c)\n");
        }
      for (i=0;i<nelray;i++)
        printf("              %6d      %25.18e            %25.18e\n",
               i,erayeq[i],trayeq[i]-tkel);
    }

  free(e);
}
/*|======================================================================|
  | SYRTHES 5.0                                       COPYRIGHT EDF 2009 |
  |======================================================================|
  | AUTEURS  : C. PENIGUEL, I. RUPP                                      |
  |======================================================================|
  | fi2teq_parall  - version parallele                                   |
  |                                                                      |
  |======================================================================| */
void fi2teq_parall(int nelray,double *tmpray,
                   struct FacForme ***listefdf,double *diagfdf,
                   double *sufray,
                   struct ProphyRay phyray,
                   double **radios,double **firay,
                   double *trayeq,double *erayeq,
                   struct PropInfini propinf,struct Clim fimpray,
                   struct Soleil soleil,
                   struct Vitre vitre,
                   struct SDparall_ray sdparall_ray,
                   double *trav0,double *trav1,double *trav2)
{
  int nb,n,m,i,j,ngfac,nel,nume,iv,jv,numbs,nf;
  double eps;
  double w1,w2,tfac;
  double ee,c2,x1,x2,x,v,v2,v4,**e,**ev,f;
  double *fdfnp1;

  fdfnp1=propinf.fdfnp1;

  if (phyray.bandespec.nb>100)
    {
      if (SYRTHES_LANG == FR)
        {
          printf("fi2teq : le nombre de bandes spectrales ne peut etre superieur a 100\n");
          printf("         dans le cas improbable ou cela pose pb contacter CP/IR \n");
        }
      else if (SYRTHES_LANG == EN)
        {
          printf("fi2teq : the number of spectral band cannot be above 100\n");
          printf("         in the very unlikely case it induces a pb contact  CP/IR \n");
        }
      syrthes_exit(1);
    }

  eps=1e-6;
  c2=1.4388e-2;
  for (i=0;i<nelray;i++) trayeq[i]=erayeq[i]=0.;

  e=(double**)malloc(phyray.bandespec.nb*sizeof(double*));
  for (n=0;n<phyray.bandespec.nb;n++) e[n]=(double*)malloc(nelray*sizeof(double));
  verif_alloue_double2d(phyray.bandespec.nb,"fi2teq",e);
  for (i=0;i<nelray;i++) for (n=0;n<phyray.bandespec.nb;n++) {e[n][i]=0;}


  /* emissivite equivalente */
  if (phyray.bandespec.nb==1)
    {
      for (i=0;i<nelray;i++) erayeq[i]=phyray.emissi[0][i];
    }
  else
    for (n=0;n<phyray.bandespec.nb;n++)
      {
        x1=c2/phyray.bandespec.borneinf[n]; x2=c2/phyray.bandespec.bornesup[n];

        for (i=0;i<nelray;i++)
          {
            v=x1/tmpray[i]; w1=wiebel(v);
            v=x2/tmpray[i]; w2=wiebel(v);
            erayeq[i] += (w2-w1)*phyray.emissi[n][i];
          }
      }


  /* temperature equivalente */
  if (phyray.bandespec.nb==1)
    {

      fdf_Y_parall(trayeq,radios[0],listefdf,diagfdf,trav0,sdparall_ray);

      for (i=0;i<nelray;i++)
        {
          if (propinf.actif)
            trayeq[i] += fdfnp1[i]*sigma*propinf.EO[0][i] * pow(propinf.TO[i]+tkel,4);

          trayeq[i]=pow(trayeq[i]/(sufray[i]*sigma),0.25);
        }
    }
  else
    {
      for (i=0;i<nelray;i++) trav2[i]=0;

      for (nb=0;nb<phyray.bandespec.nb;nb++)
        {
          fdf_Y_parall(trav1,radios[nb],listefdf,diagfdf,trav0,sdparall_ray);

          for (i=0;i<nelray;i++)
            {
              if (propinf.actif)
                {
                  x1=c2/phyray.bandespec.borneinf[nb]; x2=c2/phyray.bandespec.bornesup[nb];

                  v=x1/(propinf.TO[i]+tkel);  w1=wiebel(v);
                  v=x2/(propinf.TO[i]+tkel);  w2=wiebel(v);

                  trav1[i] += fdfnp1[i]*sigma*propinf.EO[nb][i]*pow(propinf.TO[i]+tkel,4)*(w2-w1);
                }

              trav2[i]+=trav1[i]*phyray.emissi[nb][i];
            }
        }

      for (i=0;i<nelray;i++)
        trayeq[i]=pow(trav2[i]/(sufray[i]*erayeq[i]*sigma),0.25);

    }



  /* 4.- calcul du flux de rayonnemnent (pour le post-processing) */
  for (nb=0;nb<phyray.bandespec.nb;nb++)
    {
      fdf_Y_parall(trav1,radios[nb],listefdf,diagfdf,trav0,sdparall_ray);

      for (i=0;i<nelray;i++)
        {
          if (propinf.actif)
            {
              x1=c2/phyray.bandespec.borneinf[nb]; x2=c2/phyray.bandespec.bornesup[nb];
              v=x1/(propinf.TO[i]+tkel);  w1=wiebel(v);
              v=x2/(propinf.TO[i]+tkel);  w2=wiebel(v);

              trav1[i] += fdfnp1[i]*sigma*propinf.EO[nb][i]*pow(propinf.TO[i]+tkel,4)*(w2-w1);
            }

          if (phyray.bandespec.nb==1)
            firay[nb][i]= phyray.emissi[nb][i]*(sigma*pow(tmpray[i],4)-trav1[i]/sufray[i]);
          else
            {
              v=x1/tmpray[i]; w1=wiebel(v);
              v=x2/tmpray[i]; w2=wiebel(v);
              firay[nb][i]= phyray.emissi[nb][i]*(sigma*(w2-w1)*pow(tmpray[i],4)-trav1[i]/sufray[i]);
            }
        }
    }

  /*
    if (vitre.actif)
    if (bandespec.nb==1)
    firay[0][i]= emissi[0][i]*sigma*pow(tmpray[i],4) - absorp[0][i]*ev[0][i];
    else
    for (nb=0;nb<bandespec.nb;nb++)
    for (i=0;i<nelray;i++)
    {
    v=x1/tmpray[ig]; w1=wiebel(v);
    v=x2/tmpray[ig]; w2=wiebel(v);
    firay[nb][i]= emissi[nb][i]*sigma*pow(tmpray[i],4) - absorp[n][i]*ev[n][i];
    }
    else
    if (bandespec.nb==1)
    firay[0][i]= emissi[0][i]*sigma*pow(tmpray[i],4) - absorp[0][i]*e[0][i];
    else
    for (nb=0;nb<bandespec.nb;nb++)
    for (i=0;i<nelray;i++)
    {
    v=x1/tmpray[ig]; w1=wiebel(v);
    v=x2/tmpray[ig]; w2=wiebel(v);
    firay[nb][i]= emissi[nb][i]*sigma*pow(tmpray[i],4) - absorp[nb][i]*e[nb][i];
    }
    */



  /* 5-mise a jour des facettes avec flux imposee */
  for (n=0;n<phyray.bandespec.nb;n++)
    for (i=0;i<fimpray.nelem;i++)
      {
        ngfac=fimpray.numf[i];
        firay[n][ngfac]=fimpray.val1[n][i];
      }

  /* 6-calcul de la temperature de la facette necessaire */
  for (i=0;i<fimpray.nelem;i++)
    {
      tfac=0.;
      ngfac=fimpray.numf[i];
      for (n=0;n<phyray.bandespec.nb;n++)
        tfac += (1.-phyray.emissi[n][ngfac])/phyray.emissi[n][ngfac]*fimpray.val1[n][i]
          + radios[n][i];
      tmpray[ngfac]=pow(tfac/sigma,0.25);
    }


   /* 7-calcul de la temperature de la facette necessaire pour la voute */
  double xflux;
  if (soleil.actif && (soleil.model==3 || soleil.h*180/Pi>=1.))
    {
      for (i=0;i<soleil.voute->nelem;i++)
        {
          tfac=0.;
          nf=soleil.voute->numf[i];
          for (n=0;n<phyray.bandespec.nb;n++)
            {
              xflux=soleil.voute->fluxdiffus[n];
              if (nf== soleil.voute->numf[soleil.voute->numFacePosSoleil]) xflux+=soleil.voute->fluxdirect[n];

              tfac += (1.-phyray.emissi[n][nf])/phyray.emissi[n][nf]*xflux +  radios[n][nf];
            }
          tmpray[nf]=pow(tfac/sigma,0.25);
          //printf("voute i=%d tmpray[%d]=%f radios=%f\n",i, nf,tmpray[nf],radios[0][ngfac]);

        }
    }


  if (affich.ray_fi2teq)
    {
      if (SYRTHES_LANG == FR)
        {
          printf("\n *** fi2teq rang %d : calcul des emissivites equivalentes ",syrglob_rang);
          printf(" des temperatures de rayonnement equivalentes\n");
          printf("             facette  emissivite equi      temp equivalente (degres c)\n");
        }
      else if (SYRTHES_LANG == EN)
        {
          printf("\n *** fi2teq rang %d : calculation of equivalent emissivity ",syrglob_rang);
          printf(" and equivalent radiation temperature\n");
          printf("             face  emissivity equi      equivalente temp (degres c)\n");
        }
      for (i=0;i<nelray;i++)
        printf("              %6d       %25.18e               %25.18e\n",
               i,erayeq[i],trayeq[i]-tkel);
    }

  free(e);

}
/*|======================================================================|
  | SYRTHES 5.0                                       COPYRIGHT EDF 2015 |
  |======================================================================|
  | AUTEURS  : C. PENIGUEL, I. RUPP                                      |
  |======================================================================|
  | fi2teqfdg                                                            |
  | Calcul de la temperature equivalente dans le cas des                 |
  | facteurs de Gebhart                                                  |
  |======================================================================| */
void fi2teqfdg(int nelray,double *tmpray,
               struct FacGebhart **listefdg,double **diagfdg,
               double *sufray,struct ProphyRay phyray,double **firay,
               double *trayeq,double *erayeq,
               struct PropInfini propinf,struct Clim fimpray,
               struct Soleil soleil,
               struct Vitre vitre,
               struct Fluid0d fluid0d,struct SymPer raysymper)
{
  int nb,n,m,i,j,ngfac,nel,nume,iv,jv,numbs,nc;
  double eps,ts;
  double w1,w2,tfac;
  double ee,c2,x1,x2,x,v,v2,v4,rr,**e,**ev,xx,f;
  struct FacGebhart *pfg;
  double **fdgnp1;


  fdgnp1=propinf.fdgnp1;

  if (phyray.bandespec.nb>100)
    {
      if (SYRTHES_LANG == FR)   {
          printf("fi2teqfdg : le nombre de bandes spectrales ne peut etre superieur a 100\n");
          printf("            dans le cas improbable ou cela pose pb contacter CP/IR \n");
      }
      else if (SYRTHES_LANG == EN)      {
        printf("fi2teqfdg : the number of spectral band cannot be above 100\n");
        printf("            in the very unlikely case it induces a pb contact  CP/IR \n");
      }
      syrthes_exit(1);
    }

  eps=1e-6;
  c2=1.4388e-2;

  for (i=0;i<nelray;i++) trayeq[i]=erayeq[i]=0.;



  if (phyray.bandespec.nb==1)
    {
      for (i=0;i<nelray;i++) erayeq[i]=phyray.emissi[0][i]*(1-diagfdg[0][i]);
    }
  else
    for (n=0;n<phyray.bandespec.nb;n++)
      {
        x1=c2/phyray.bandespec.borneinf[n]; x2=c2/phyray.bandespec.bornesup[n];

        for (i=0;i<nelray;i++)
          {
            v=x1/tmpray[i]; w1=wiebel(v);
            v=x2/tmpray[i]; w2=wiebel(v);
            erayeq[i] += (w2-w1)*phyray.emissi[n][i]*(1-diagfdg[n][i]);
          }
      }



  /* --- calcul de la temperature equivalente  ----- */
  /* ----------------------------------------------- */
  if (phyray.bandespec.nb==1)
    {
      for (i=0;i<nelray;i++)
        {
          rr=0;
          pfg=listefdg[i];
          for (m=0;m<listefdg[i]->nb;m++)  rr+=pfg->fdg[0][m]*pow(tmpray[pfg->numenface[m]],4);

          //      if (propinf.actif) rr += fdgnp1[0][i]*xj[0];
          if (lfluid0d==2 && fluid0d.caviteF[i]>-1)  rr += fluid0d.fdgv[0][i]*pow(fluid0d.cavite[fluid0d.caviteF[i]].Tf+tkel,4);
          rr /= (1-diagfdg[0][i]);
          trayeq[i]=pow(rr,0.25);
        }
    }
  else
    {
      for (i=0;i<nelray;i++)
        {
          rr=0.;
          pfg=listefdg[i];
          for (nb=0;nb<phyray.bandespec.nb;nb++)
            {
              x1=c2/phyray.bandespec.borneinf[nb]; x2=c2/phyray.bandespec.bornesup[nb];
              xx=0;

              for (m=0;m<listefdg[i]->nb;m++)
                {
                  xx+=pfg->fdg[nb][m]*pow(tmpray[pfg->numenface[m]],4);
                  v=x1/tmpray[pfg->numenface[m]]; w1=wiebel(v);
                  v=x2/tmpray[pfg->numenface[m]]; w2=wiebel(v);
                  xx *= (w2-w1);
                }

              //              if (propinf.actif) xx += fdgnp1[i]*xj[nb];
              if (lfluid0d==2 && fluid0d.caviteF[i]>-1)  xx += fluid0d.fdgv[nb][i]*pow(fluid0d.cavite[fluid0d.caviteF[i]].Tf+tkel,4);

              rr += xx;
            }
          trayeq[i]=pow(rr/erayeq[i],0.25);
        }
    }

  /* modele fluide 0D gaz gris par bande : calcul du terme source */
  if (lfluid0d==2 && lray_fdf==2)
    {
      for (n=0;n<fluid0d.nb_cavite;n++) fluid0d.cavite[n].sourceray=0;

      if (phyray.bandespec.nb==1)
        for (i=0;i<nelray;i++)
          {
            nc=fluid0d.caviteF[i];
            fluid0d.cavite[nc].sourceray += phyray.emissi[0][i]*sufray[i]*sigma*(pow(tmpray[i],4)-pow(fluid0d.cavite[nc].Tf+tkel,4))*fluid0d.fdgv[0][i];
          }
      else
        for (i=0;i<nelray;i++)
          {
            nc=fluid0d.caviteF[i];
            for (nb=0;nb<phyray.bandespec.nb;nb++)
              {
                x1=c2/phyray.bandespec.borneinf[nb]; x2=c2/phyray.bandespec.bornesup[nb];
                v=x1/tmpray[i]; w1=wiebel(v);
                v=x2/tmpray[i]; w2=wiebel(v);
                fluid0d.cavite[nc].sourceray += (phyray.emissi[nb][i]*sufray[i]*sigma*(pow(tmpray[i],4)-pow(fluid0d.cavite[nc].Tf+tkel,4))*fluid0d.fdgv[nb][i])*(w2-w1);
              }
          }

      if (raysymper.nbsym || raysymper.nbper)
        for (n=0;n<fluid0d.nb_cavite;n++)
          {
            ts=fluid0d.cavite[n].sourceray;
            if (raysymper.nbsym) fluid0d.cavite[nc].sourceray = ts*pow(2,raysymper.nbsym);
            if (raysymper.nbper) fluid0d.cavite[nc].sourceray += ts*raysymper.nbper;
          }
    }



  /* 4.- calcul du flux net de rayonnemnent (pour le post-processing) */
  if (phyray.bandespec.nb==1)
    {
      nb = 0;
      for (i=0;i<nelray;i++) firay[0][i]= phyray.emissi[nb][i]*sigma*sufray[i] * (pow(tmpray[i],4)-pow(trayeq[i],4));
    }
  else
    {
      for (nb=0;nb<phyray.bandespec.nb;nb++)
        for (i=0;i<nelray;i++)
          {
            rr=0;
            for (m=0;m<listefdg[i]->nb;m++)
              {
                rr+=listefdg[i]->fdg[nb][m]*pow(tmpray[listefdg[i]->numenface[m]],4);
                v=x1/tmpray[pfg->numenface[m]]; w1=wiebel(v);
                v=x2/tmpray[pfg->numenface[m]]; w2=wiebel(v);
                rr *= (w2-w1);
            }
            v=x1/tmpray[i]; w1=wiebel(v);
            v=x2/tmpray[i]; w2=wiebel(v);
            rr+=diagfdg[nb][i]*pow(tmpray[i],4) * (w2-w1);

            firay[nb][i]= phyray.emissi[nb][i]*sufray[i]*(w2-w1)*pow(tmpray[i],4) - rr;
          }
    }



  /* 5-mise a jour des facettes avec flux imposee */
  for (n=0;n<phyray.bandespec.nb;n++)
    for (i=0;i<fimpray.nelem;i++)
      {
        ngfac=fimpray.numf[i];
        firay[n][ngfac]=fimpray.val1[n][i];
      }

  /* 6-calcul de la temperature de la facette necessaire */

  // pour l'instant on ne sait pas faire a partir des fdg ?????

  //  for (i=0;i<fimpray.nelem;i++)
  // {
  //  tfac=0.;
  //  ngfac=fimpray.numf[i];
  //  for (n=0;n<phyray.bandespec.nb;n++)
  //    tfac += (1.-phyray.emissi[n][ngfac])/phyray.emissi[n][ngfac]*fimpray.val1[n][i]
  //      + radios[n][i];
  //  tmpray[ngfac]=pow(tfac/sigma,0.25);
  //}


  if (affich.ray_fi2teq)
    {
      if (SYRTHES_LANG == FR)
        {
          printf("\n *** fi2teqfdg : calcul des emissivites et");
          printf(" des temperatures de rayonnement equivalentes\n");
          printf("                facette    emissivite equi   temp equi (degres c)\n");
        }
      else if (SYRTHES_LANG == EN)
        {
          printf("\n *** fi2teqfdg : calculation of equivalent emissivity ");
          printf(" and equivalent radiation temperature\n");
          printf("                face    emissivity equi   temp equi (degres c)\n");
        }
      for (i=0;i<nelray;i++)
        printf("              %6d      %25.18e            %25.18e\n",
               i,erayeq[i],trayeq[i]-tkel);
    }

}
/*|======================================================================|
  | SYRTHES 5.0                                       COPYRIGHT EDF 2015 |
  |======================================================================|
  | AUTEURS  : C. PENIGUEL, I. RUPP                                      |
  |======================================================================|
  | fi2teqfdg_parall  - version parallele de fi2teqfdg                   |
  |                                                                      |
  |======================================================================| */
void fi2teqfdg_parall(int nelray,double *tmpray,
                      struct FacGebhart ***listefdg,double **diagfdg,
                      double *sufray,
                      struct ProphyRay phyray,
                      double **firay,
                      double *trayeq,double *erayeq,
                      struct PropInfini propinf,struct Clim fimpray,
                      struct Soleil soleil,
                      struct Vitre vitre,struct Fluid0d fluid0d,struct SymPer raysymper,
                      struct Travail trav,struct SDparall_ray sdparall_ray)

{
  int nb,n,m,i,j,ngfac,nel,nume,iv,jv,numbs,nc;
  double eps,ts;
  double w1,w2,tfac;
  double xx,ee,c2,x1,x2,x,v,v2,v4,**e,**ev,f;
  double *fdfnp1;
  double *trav0,*trav1,*trav2,*trav3,*trav4;


  trav0=trav.tab[0];
  trav1=trav.tab[1];
  trav2=trav.tab[2];
  trav3=trav.tab[3];
  trav4=trav.tab[4];

#ifdef _SYRTHES_MPI_

  fdfnp1=propinf.fdfnp1;

  if (phyray.bandespec.nb>100)
    {
      if (SYRTHES_LANG == FR)
        {
          printf("fi2teqfdg_parall : le nombre de bandes spectrales ne peut etre superieur a 100\n");
          printf("                   dans le cas improbable ou cela pose pb contacter CP/IR \n");
        }
      else if (SYRTHES_LANG == EN)
        {
          printf("fi2teqfdg_parall : the number of spectral band cannot be above 100\n");
          printf("                   in the very unlikely case it induces a pb contact  CP/IR \n");
        }
      syrthes_exit(1);
    }

  eps=1e-6;
  c2=1.4388e-2;
  for (i=0;i<nelray;i++) trayeq[i]=erayeq[i]=0.;



  if (phyray.bandespec.nb==1)
    {
      for (i=0;i<nelray;i++) erayeq[i]=phyray.emissi[0][i]*(1-diagfdg[0][i]);
    }
  else
    for (n=0;n<phyray.bandespec.nb;n++)
      {
        x1=c2/phyray.bandespec.borneinf[n]; x2=c2/phyray.bandespec.bornesup[n];

        for (i=0;i<nelray;i++)
          {
            v=x1/tmpray[i]; w1=wiebel(v);
            v=x2/tmpray[i]; w2=wiebel(v);
            erayeq[i] += (w2-w1)*phyray.emissi[n][i]*(1-diagfdg[n][i]);
          }
      }

  /* --- calcul de la temperature equivalente  ----- */
  /* ----------------------------------------------- */

  for (i=0;i<nelray;i++) trav1[i]=pow(tmpray[i],4);


  if (phyray.bandespec.nb==1)
    {
      fdg_Y_parall(trayeq,trav1,listefdg,diagfdg,0,trav0,sdparall_ray);

      for (i=0;i<nelray;i++)
        {
          //      if (propinf.actif) trayeq[i] += fdfnp1[i]*xj[0];
          if (lfluid0d==2 && fluid0d.caviteF[i]>-1)  trayeq[i] += fluid0d.fdgv[0][i]*pow(fluid0d.cavite[fluid0d.caviteF[i]].Tf+tkel,4);
          trayeq[i]/=(1-diagfdg[0][i]);

          trayeq[i]=pow(trayeq[i],0.25);
        }
    }
  else
    {
      for (i=0;i<nelray;i++) trav3[i]=0;

      for (nb=0;nb<phyray.bandespec.nb;nb++)
        {
          x1=c2/phyray.bandespec.borneinf[nb]; x2=c2/phyray.bandespec.bornesup[nb];
          xx=0;

          fdg_Y_parall(trav2,trav1,listefdg,diagfdg,nb,trav0,sdparall_ray);

          for (i=0;i<nelray;i++)
            {
              v=x1/tmpray[i]; w1=wiebel(v);
              v=x2/tmpray[i]; w2=wiebel(v);
              trav3[i] += trav2[i]*(w2-w1);

              //              if (propinf.actif) trav1[i] += fdfnp1[i]*xj[nb];
              if (lfluid0d==2 && fluid0d.caviteF[i]>-1)  trav3[i] += fluid0d.fdgv[0][i]*pow(fluid0d.cavite[fluid0d.caviteF[i]].Tf+tkel,4);
            }

        }


      for (i=0;i<nelray;i++) trayeq[i]=pow(trav3[i]/erayeq[i],0.25);

    }

  /* modele fluide 0D gaz gris par bande : calcul du terme source */
  if (lfluid0d==2 && lray_fdf==2)
    {
      for (n=0;n<fluid0d.nb_cavite;n++) fluid0d.cavite[n].sourceray=trav3[n]=0;

      if (phyray.bandespec.nb==1)
        for (i=0;i<nelray;i++)
          {
            nc=fluid0d.caviteF[i];
            trav3[nc] += phyray.emissi[0][i]*sufray[i]*sigma*(pow(tmpray[i],4)-pow(fluid0d.cavite[nc].Tf+tkel,4))*fluid0d.fdgv[0][i];
          }
      else
        for (i=0;i<nelray;i++)
          {
            nc=fluid0d.caviteF[i];
            for (nb=0;nb<phyray.bandespec.nb;nb++)
              {
                x1=c2/phyray.bandespec.borneinf[nb]; x2=c2/phyray.bandespec.bornesup[nb];
                v=x1/tmpray[i]; w1=wiebel(v);
                v=x2/tmpray[i]; w2=wiebel(v);
                trav3[nc] += (phyray.emissi[nb][i]*sufray[i]*sigma*(pow(tmpray[i],4)-pow(fluid0d.cavite[nc].Tf+tkel,4))*fluid0d.fdgv[nb][i])*(w2-w1);
              }
          }

      for (n=0;n<fluid0d.nb_cavite;n++) trav4[n]=0;
      MPI_Allreduce(trav3,trav4,fluid0d.nb_cavite,MPI_DOUBLE,MPI_SUM,sdparall_ray.syrthes_comm_ray);
      for (n=0;n<fluid0d.nb_cavite;n++) fluid0d.cavite[n].sourceray=trav4[n];

      if (raysymper.nbsym || raysymper.nbper)
        for (n=0;n<fluid0d.nb_cavite;n++)
          {
            ts=fluid0d.cavite[n].sourceray;
            if (raysymper.nbsym) fluid0d.cavite[n].sourceray = ts*pow(2,raysymper.nbsym);
            if (raysymper.nbper) fluid0d.cavite[n].sourceray += ts*raysymper.nbper;
          }
    }



  /* 4.- calcul du flux de rayonnemnent (pour le post-processing) */
  for (nb=0;nb<phyray.bandespec.nb;nb++)
    {
      fdg_Y_parall(trav2,trav1,listefdg,diagfdg,nb,trav0,sdparall_ray);

      for (i=0;i<nelray;i++)
        {
          //      if (propinf.actif) trav1[i] += fdfnp1[i]*xj[nb];
          if (phyray.bandespec.nb==1)
            firay[nb][i]= phyray.emissi[nb][i]*(sigma*pow(tmpray[i],4)-trav2[i]/sufray[i]);
          else
            {
              v=x1/tmpray[i]; w1=wiebel(v);
              v=x2/tmpray[i]; w2=wiebel(v);
              firay[nb][i]= phyray.emissi[nb][i]*(sigma*(w2-w1)*pow(tmpray[i],4)-trav2[i]/sufray[i]);
            }
        }
    }


  /* 5-mise a jour des facettes avec flux imposee */
  for (n=0;n<phyray.bandespec.nb;n++)
    for (i=0;i<fimpray.nelem;i++)
      {
        ngfac=fimpray.numf[i];
        firay[n][ngfac]=fimpray.val1[n][i];
      }

  /* 6-calcul de la temperature de la facette necessaire */
 // pour l'instant on ne sait pas faire a partir des fdg ?????
  /* for (i=0;i<fimpray.nelem;i++) */
  /*   { */
  /*     tfac=0.; */
  /*     ngfac=fimpray.numf[i]; */
  /*     for (n=0;n<phyray.bandespec.nb;n++)  */
  /*    tfac += (1.-phyray.emissi[n][ngfac])/phyray.emissi[n][ngfac]*fimpray.val1[n][i]  */
  /*      + radios[n][i]; */
  /*     tmpray[ngfac]=pow(tfac/sigma,0.25); */
  /*   } */


  if (affich.ray_fi2teq)
    {
      if (SYRTHES_LANG == FR)
        {
          printf("\n *** fi2teq rang %d : calcul des emissivites equivalentes ",syrglob_rang);
          printf(" des temperatures de rayonnement equivalentes\n");
          printf("             facette  emissivite equi      temp equivalente (degres c)\n");
        }
      else if (SYRTHES_LANG == EN)
        {
          printf("\n *** fi2teq rang %d : calculation of equivalent emissivity ",syrglob_rang);
          printf(" and equivalent radiation temperature\n");
          printf("             face  emissivity equi      equivalente temp (degres c)\n");
        }
      for (i=0;i<nelray;i++)
        printf("              %6d       %25.18e               %25.18e\n",
               i,erayeq[i],trayeq[i]-tkel);
    }

  /* penser a faire un free pour e et ev ?????????????????????? */

#endif
}
/*|======================================================================|
  | SYRTHES 5.0                                       COPYRIGHT EDF 2009 |
  |======================================================================|
  | AUTEURS  : C. PENIGUEL, I. RUPP                                      |
  |======================================================================|
  | gausei                                                               |
  |                                                                      |
  |======================================================================| */
void rrayrc_parall(int nelray,struct FacForme ***listefdf,double *diagfdf,
                   double *sufray,
                   double *rho,double *epr,
                   struct MatriceRay matrray,int *niter,
                   double *trav0,struct SDparall_ray sdparall_ray)
{
  int n,ngfac;
  double aa,bb,cc,d,ee,alfa,aalfa,tau,alfam1,alftau;
  double aux1,aux2,aux,denom;
  int i,j,nitsmo;
  double x0,resnor,prsca1;
  double epsis,zero,epssmo;
  struct FacForme *pff;
  double *x,*b,*xm1,*gd,*res,*z,*di,*resm1;


  x     = matrray.x;
  b     = matrray.b;
  xm1   = matrray.xm1;
  gd    = matrray.gd;
  res   = matrray.res;
  z     = matrray.z;
  di    = matrray.di;
  resm1 = matrray.resm1;


  zero=0.;
  nitsmo=100;
  epssmo=1.e-12;
  n=0;

  /* 1- initialisation des vecteurs auxiliaires */
  for (i=0;i<nelray;i++)
    {
      xm1[i]=x[i]=epr[i];
      b[i]=epr[i];
      di[i]=sufray[i]-rho[i]*diagfdf[i];
    }

  /* norme du second membre */
  prosca_ray_parall(&prsca1,x,x,sdparall_ray);
  x0=sqrt(prsca1);

  if (x0< 1.e-20 ) x0=1.e-6;
  epsis=1.e-4 * x0;


  fdf_Y_parall(res,x,listefdf,diagfdf,trav0,sdparall_ray);


  for (i=0;i<nelray;i++)
    res[i] = -rho[i]*res[i] + di[i]*x[i] - b[i] + rho[i]*x[i]*diagfdf[i];

  prosca_ray_parall(&prsca1,res,res,sdparall_ray);
  resnor=sqrt(prsca1);
  /*   printf(" >>>> ray=resnor_ini par : %25.18e \n",resnor); */



  if (sdparall_ray.rangray==0)
    {
      if (SYRTHES_LANG == FR)
        printf(" *** RRAYRC : ITERATION   PRECISION RELATIVE  PRECISION ABSOLUE\n");
      else if (SYRTHES_LANG == EN)
        printf(" *** RRAYRC : ITERATION   RELATIVE PRECISION  ABSOLUTE PRECISION\n");
    }

  if ( resnor<=epsis && resnor<=epsgcs_rad*sqrt((double)(sdparall_ray.nelraytot)) )
    {
      if (sdparall_ray.rangray==0)
        printf("                 %4d         %12.5e                %12.5e\n",
               n,resnor/x0,resnor/sqrt((double)(sdparall_ray.nelraytot)));
      *niter=n;
    }

  else
    {
      for (i=0;i<nelray;i++) resm1[i]=res[i];
      n=0;

      do
        {
          n++;
          for (i=0;i<nelray;i++) gd[i]=res[i]/di[i];

          fdf_Y_parall(z,gd,listefdf,diagfdf,trav0,sdparall_ray);

          for (i=0;i<nelray;i++)
            {z[i] = -rho[i]*z[i] + di[i]*gd[i] + rho[i]*gd[i]*diagfdf[i];    trav0[i]=res[i]-resm1[i];}

          aa=0.;bb=0.;cc=0.;d =0.;ee=0.;
          prosca_ray_parall(&aa,z,res,sdparall_ray);
          prosca_ray_parall(&bb,z,resm1,sdparall_ray);
          prosca_ray_parall(&cc,res,trav0,sdparall_ray);
          prosca_ray_parall(&d,resm1,trav0,sdparall_ray);
          prosca_ray_parall(&ee,z,z,sdparall_ray);

          denom=(cc-d )*ee-(aa-bb)*(aa-bb);
          if(fabs(denom)< 1.e-20) alfa=1.;
          else alfa=((aa-bb)*bb-d *ee)/denom;
          aalfa=fabs(alfa);
          if (aalfa<1e-20 || fabs(aalfa-1.)< 1.e-20) {alfa=1.; tau=aa/ee;}
          else tau=aa/ee + (1.-alfa)/alfa * bb/ee;
          alfam1 =1.-alfa;
          alftau =-alfa*tau;
          for (i=0;i<nelray;i++)
            {
              aux1    =res[i];
              aux2    =x[i];
              res[i]  =alfa*aux1+alfam1*resm1[i]+alftau*z[i];
              resm1[i]=aux1;
              x[i]    =alfa*aux2+alfam1*xm1[i]+alftau*gd[i];
              xm1[i]  =aux2;
            }

          prosca_ray_parall(&prsca1,res,res,sdparall_ray);
          resnor=sqrt(prsca1);

/*        if (sdparall_ray.rangray==0) printf("rrayrc_parall : it=%d   resnor=%25.18e \n",n,resnor); */

          if (sdparall_ray.rangray==0 && n%10==0)
            printf("                 %4d         %12.5e        %12.5e\n",
                   n,resnor/x0,resnor/sqrt((double)(sdparall_ray.nelraytot)));

        }
      while(!( (resnor<=epsis &&  resnor<=epsgcs_rad*sqrt((double)(sdparall_ray.nelraytot))) || n>=nitsmo));


      if (sdparall_ray.rangray==0 && n%10!=0)
        printf("                 %4d         %12.5e        %12.5e\n",
               n,resnor/x0,resnor/sqrt((double)(sdparall_ray.nelraytot)));
      *niter=n;
    }

}
/*|======================================================================|
  | SYRTHES 5.0                                       COPYRIGHT EDF 2009 |
  |======================================================================|
  | AUTEURS  : I. RUPP, C. PENIGUEL                                      |
  |======================================================================|
  | fdf_Y_parall                                                         |
  |         Calcul matrice vecteur   FDF.J en parallele (rayonnement)    |
  |        RES = fdf . Y   (AVEC la diagonale)                           |
  |======================================================================|*/
void fdf_Y_parall(double *res,double *y,
                  struct FacForme ***listefdf,double *diagfdf,
                  double *trav0,struct SDparall_ray sdparall_ray )
{
#ifdef _SYRTHES_MPI_
  int i,n,nn,ideb,ifin,nelemloc,indic;
  struct FacForme *pff;

  ideb=sdparall_ray.ieledeb[sdparall_ray.rangray];
  ifin=sdparall_ray.ieledeb[sdparall_ray.rangray+1];
  nelemloc=sdparall_ray.nelrayloc[sdparall_ray.rangray];



  /* multiplication en local */
  /* ----------------------- */
  for (i=0;i<nelemloc;i++) res[i]=0.;

  for (i=0;i<nelemloc;i++)
    for (nn=0;nn<listefdf[sdparall_ray.rangray][i]->nb;nn++)
      res[i] += listefdf[sdparall_ray.rangray][i]->fdf[nn] * y[listefdf[sdparall_ray.rangray][i]->numenface[nn]];


  /* on complete avec ce qui provient des autres processeurs */
  /* ------------------------------------------------------- */

  /* Pour les fdf en raytracing c'est particulier car la matrice n'est pas forcement symetrique                 */
  /* le proc i peut avoir des fdf avec le proc j ALORS QUE le proc j n'en a pas avec le proc i !  */
    for (n=0;n<sdparall_ray.npartsray;n++)
    {
      if (n!=sdparall_ray.rangray)
        {
          indic = sdparall_ray.fdfproc[n];

          /* est-ce-que le proc n voit des facettes de mon proc courant ?  ==> recuperer de fdfproc du proc pour mon rang */
          /* j'envoie l'indicateur que je vois des facettes du proc n                                                     */
          /* si sdparall_ray.fdfproc[n] ==> je vois des facettes du proc n ==> il y a des facettes a recuperer du proc n  */

          MPI_Sendrecv_replace(&indic,1,MPI_INT,n,0,n,0,sdparall_ray.syrthes_comm_ray,&status);


          /* si l'echange est symetrique, on echange les donnees des procs */
          if (sdparall_ray.fdfproc[n] && indic==1)
            {
              for (i=0;i<nelemloc;i++) trav0[i]=y[i];
              for (i=nelemloc;i<sdparall_ray.nelrayloc_max;i++) trav0[i]=0.;

              MPI_Sendrecv_replace(trav0,sdparall_ray.nelrayloc_max,MPI_DOUBLE,n,0,n,0,sdparall_ray.syrthes_comm_ray,&status);

              /* on fait la multiplication avec ce que l'on a recu */
              for (i=0;i<nelemloc;i++)
                for (nn=0;nn<listefdf[n][i]->nb;nn++)
                  res[i] += listefdf[n][i]->fdf[nn] * trav0[listefdf[n][i]->numenface[nn]];
            }

          /* si indic=1 ==> le proc n a besoin de mon Y local */
          else if (indic==1)
            {
              for (i=0;i<nelemloc;i++) trav0[i]=y[i];
              for (i=nelemloc;i<sdparall_ray.nelrayloc_max;i++) trav0[i]=0.;
              MPI_Send(trav0,sdparall_ray.nelrayloc_max,MPI_DOUBLE,n,0,sdparall_ray.syrthes_comm_ray);
            }


          else if (sdparall_ray.fdfproc[n]) /* j'ai besoin du Y du proc n */
            {
              MPI_Recv(trav0,sdparall_ray.nelrayloc_max,MPI_DOUBLE,n,0,sdparall_ray.syrthes_comm_ray,&status);

              /* on fait la multiplication avec ce que l'on a recu */
              for (i=0;i<nelemloc;i++)
                for (nn=0;nn<listefdf[n][i]->nb;nn++)
                  res[i] += listefdf[n][i]->fdf[nn] * trav0[listefdf[n][i]->numenface[nn]];
            }
        }
    }

  /* ajout de la diagonale */
  for (i=0;i<nelemloc;i++) res[i] += diagfdf[i]*y[i];

#endif
}