*** /home/muntwiler_m/mnt/pearl_data/software/edac/edac_all.cpp	2011-04-14 23:38:44.000000000 +0200
--- edac_all.cpp	2018-02-05 17:30:17.347373088 +0100
***************
*** 3085,3090 ****
--- 3085,3091 ----
  numero Vxc_Barth_Hedin(numero den, numero denup)
  {
    if(den<1e-10)  return 0;
+   if(denup<0)  denup=0;
    numero rs=1/pow(4*pi/3*den, 1/3.0);
    numero x=denup/den;
    numero alpha0=pow(4/(9*pi), 1/3.0);
***************
*** 10117,10122 ****
--- 10118,10124 ----
    void    scan_imfp(char *name);      
    void    scan_imfp(FILE *fout);      
    numero  iimfp_TPP(numero kr);
+   numero  iimfp_SD(numero kr);
    numero  TPP_rho, TPP_Nv, TPP_Ep, TPP_Eg;
    numero  screening_length;   
    int     scattering_so;      
***************
*** 10230,10235 ****
--- 10232,10238 ----
    
    int     n_th;
    int     n_fi;
+   int     n_ang;
    numero  *th, *fi;
    
    numero  *th_out,          
***************
*** 10239,10244 ****
--- 10242,10248 ----
    void free(void);
    void init_th(numero thi, numero thf, int nth);
    void init_phi(numero fii, numero fif, int nfi);
+   void read_angles(FILE *fin, char *my_file);
    void init_refraction(
              numero refraction);  
    void init_transmission(        
***************
*** 10905,10942 ****
  }
  numero calculation::IIIthfi(int no, numero th, numero fi)
  {
!   numero ii,ii0, xth,xfi;
!   numero th0=final.th[0], th1=final.th[final.n_th-1];
    numero fi0=final.fi[0], fi1=final.fi[final.n_fi-1];
    int ith,ifi,ij;
    while(fi<0)    fi+=2*pi;
    while(fi>2*pi) fi-=2*pi;
!   xth=(final.n_th-1)*(th-th0)/(th1-th0);   ith=int(floor(xth-0.001));
!   xfi=(final.n_fi-1)*(fi-fi0)/(fi1-fi0);   ifi=int(floor(xfi-0.001));
    if(ifi==-1)  ifi=0;
    if(0<=ith && ith<final.n_th-1 && 0<=ifi && ifi<final.n_fi-1) {
      ij=no*n_ang+ith*final.n_fi+ifi;
!     ii0=III0[ij];
!     ii=ii0 + (xth-ith)*(III0[ij+final.n_fi]-ii0) + (xfi-ifi)*(III0[ij+1]-ii0);
    }  else  ii=0;
    if(ii<0)  ii=0;
    return ii;
  }
! numero calculation::IIIave(int no, numero th, numero fi)
  {
!   if(thave<=1e-6)  return IIIthfi(no,th,fi);
    int  i,j, nn=10, mm=50;
!   numero  tth,ffi,val=0, r,f;
    for(i=0; i<nn; i++)
    for(j=0; j<mm; j++)  {
!     r=i*thave/nn;
!     f=j*2*pi/mm;
!     tth=th+r*cos(f);
!     ffi=fi+r*sin(f)/cos(pi*th/180);
!     val+=IIIthfi(no,tth,ffi);
    }
!   return val/(nn*mm);
  }
  void calculation::write_ang(FILE *fout_, int ik)          
  {
    int no,nno,i,j;        
--- 10909,10963 ----
  }
  numero calculation::IIIthfi(int no, numero th, numero fi)
  {
!   numero ii,xth,xfi;
!   numero th0=final.th_out[0], th1=final.th_out[final.n_th-1];
    numero fi0=final.fi[0], fi1=final.fi[final.n_fi-1];
    int ith,ifi,ij;
    while(fi<0)    fi+=2*pi;
    while(fi>2*pi) fi-=2*pi;
!   xth=(final.n_th-1)*(th-th0)/(th1-th0);   ith=int(floor(xth-0.0001));
!   xfi=(final.n_fi-1)*(fi-fi0)/(fi1-fi0);   ifi=int(floor(xfi-0.0001));
    if(ifi==-1)  ifi=0;
+   if(ith==-1)  ith=0;
    if(0<=ith && ith<final.n_th-1 && 0<=ifi && ifi<final.n_fi-1) {
      ij=no*n_ang+ith*final.n_fi+ifi;
!     xth = xth-ith;
!     xfi = xfi-ifi;
!     ii=III0[ij]*(1-xth)*(1-xfi) + III0[ij+final.n_fi]*xth*(1-xfi) + III0[ij+1]*(1-xth)*xfi + III0[ij+final.n_fi+1]*xth*xfi;
    }  else  ii=0;
    if(ii<0)  ii=0;
    return ii;
  }
! numero calculation::IIIave(int no, numero th, numero ph)
  {
!   if(thave<=1e-6)  return IIIthfi(no,th,ph);
    int  i,j, nn=10, mm=50;
!   numero  tth, ffi, val=0, th1, ph1, cf, nw=0, x0, y0, z0, x1, y1, th2, ph2;
    for(i=0; i<nn; i++)
    for(j=0; j<mm; j++)  {
!     th1=i*2*thave/(nn-1);            //2*sigma range
!     ph1=j*2*pi/mm;
!     //rotation of (001) around Y by th1 and around Z by ph1
!     x0 = cos(ph1)*sin(th1);
!     y1 = sin(ph1)*sin(th1);
!     z0 = cos(th1);
!     //rotation around Y by th
!     x1 = x0*cos(th) + z0*sin(th);
!     z0 = -x0*sin(th) + z0*cos(th);
!     //rotation around Z by ph
!     x0 = x1*cos(ph) - y1*sin(ph);
!     y0 = x1*sin(ph) + y1*cos(ph);
! 
!     th2 = acos(z0);
!     ph2 = atan2(y0,x0);
!    
!     cf = exp(-(th1*th1/thave*thave))*(i>0.1?i:0.1);        //gauss weight * radial weight
!     nw += cf;                                    //sum of weights
!     val+=IIIthfi(no,th2,ph2)*cf;
    }
!   return val/nw;
  }
+ 
  void calculation::write_ang(FILE *fout_, int ik)        
  {
    int no,nno,i,j;      
***************
*** 10961,10967 ****
        for(no=0; no<nno; no++)          
        for(i=0; i<final.n_th; i++)
        for(j=0; j<final.n_fi; j++)
!         III[no*n_ang+i*final.n_fi+j]=IIIave(no, final.th[i], final.fi[j]);
        delete [] III0;
      }
      for(i=0; i<final.n_th; i++)
--- 10982,10988 ----
        for(no=0; no<nno; no++)        
        for(i=0; i<final.n_th; i++)
        for(j=0; j<final.n_fi; j++)
!         III[no*n_ang+i*final.n_fi+j]=IIIave(no, final.th_out[i], final.fi[j]);
        delete [] III0;
      }
      for(i=0; i<final.n_th; i++)
***************
*** 12485,12490 ****
--- 12506,12512 ----
      else {
        kr=sqrt(sqr(calc.k[ik])+2*V0);            
        if(iimfp_flag==0)  ki=iimfp.val(kr)/2;    
+       else if(iimfp_flag==3) ki=iimfp_SD(kr)/2;
        else               ki=iimfp_TPP(kr)/2;    
        set_k(complex(kr,ki));
    } }  else  if(calc.k_flag==2)  set_k(calc.kc[ik]);
***************
*** 12507,12512 ****
--- 12529,12540 ----
    numero imfp=E/(TPP_Ep*TPP_Ep*(beta*log(gamma*E)-C/E+D/(E*E)))/a0_au;
    return  1/imfp;
  }
+ numero propagation::iimfp_SD(numero kr)
+ {
+   numero E=sqr(kr)/2*au_eV;
+   numero imfp = (1.43e3/sqr(E) + 0.54*sqrt(E))/a0_au;
+   return  1/imfp;
+ }
  void propagation::scan_imfp(char *name)
  {
    FILE *fout=NULL;                           
***************
*** 13202,13208 ****
  }
  final_state::final_state(void)
  {
!   n_th=n_fi=0;
    n_1=n_2=0;
    Ylm0_th_flag=Ylm0_fi_flag=0;
    mesh_flag=0;
--- 13230,13236 ----
  }
  final_state::final_state(void)
  {
!   n_th=n_fi=n_ang=0;
    n_1=n_2=0;
    Ylm0_th_flag=Ylm0_fi_flag=0;
    mesh_flag=0;
***************
*** 13233,13238 ****
--- 13261,13289 ----
      if(n_fi==1)  fi[0]=fii;
      else  for(j=0; j<n_fi; j++)  fi[j]=fii+j*(fif-fii)/(n_fi-1);
  } }
+ void final_state::read_angles(FILE *fin, char *my_file)
+ {
+   FILE *fang;  int i, nang;
+   if(!strcmpC(my_file,"inline"))  fang=fin;
+   else                            fang=open_file(foutput,my_file,"r");
+   nang=read_int(fang);
+   free_mesh();
+   if(nang>1) {
+     delete [] th;  delete [] th_out;  delete [] transmission; delete [] fi;
+     n_th=nang;
+     th=new numero [n_th];                 
+     th_out=new numero [n_th];             
+     transmission=new numero [n_th];       
+     n_fi=nang;
+     fi=new numero [n_fi];                   
+     for(i=0; i<nang; i++) {
+       th[i]=th_out[i]=read_numero(fang);
+       transmission[i]=1;
+       fi[i]=read_numero(fang);
+     }
+   }
+   if(strcmpC(my_file,"inline"))  fclose(fang);
+ }
  void final_state::init_refraction(numero refraction)
  {
    int i;
***************
*** 14743,14748 ****
--- 14794,14801 ----
                             || scat.TPP_Ep<=0 || scat.TPP_Eg<0)
            on_error(foutput,"(input) imfp TPP-2M", "wrong parameters");
          scat.iimfp_flag=1;
+       } else if(!strcmpC(name,"SD-UC")) {
+     scat.iimfp_flag=3;
        } else {
          scat.read_imfp(fprog,name);
          scat.iimfp_flag=0;
***************
*** 15162,15164 ****
--- 15215,15224 ----
    fprintf(foutput,"That's all, folks!\n");
    return 0;
  }
+ int run_script(char *scriptfile)
+ {
+   particle_type=electrones;
+   init_fact();
+   electron.program(scriptfile);
+   return 0;
+ }