/* Written by

Eugeny B. Shybanov,
Department of Optics,
Marine Hydrophysical Institute,
Kapitanskaya 2, Sevastopol,
Crimea, Ukraine, 99000
e-mail: <lee@alpha.mhi.iuf.net>

This translation from Pascal possibly needs some debugging if you change some constants, but it is OK if you don't touch it.

This code is based on:

*/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
/*
#include <process.h>
*/
// #include <timer.cpp>
#include <math.h>
#define Pi 3.14159265358979312
#define pir180 1.74532925199432955e-2
#define ABSize 1024
#define NReper 1024
#define kl 7.0
#define NFiMax 1801
#define HalfFiMax (NFiMax+1)>>1

typedef double AngleType[NFiMax];
typedef double HalfAngleType[HalfFiMax];
typedef double PTType[HalfFiMax*2];
typedef double CType[2];
typedef CType QwType[ABSize];
typedef CType ReperType[NReper+1];
typedef double RBType8[ABSize];
typedef double fun(double x);

HalfAngleType PiM1;
HalfAngleType Pi0;
HalfAngleType mu;
int flagf=0;
char *sdf;
double a1=0,a2=0,a3=0;
int nfi,show;
long int np;
double r1,r2,la,mr,mi,n0;
FILE *otx;

//----------------------------------------------------------
double proba(double x,fun yx)
{

return yx(x);
};
//----------------------------------------
double Yunge(double x)
{
char s[20];
if (!flagf)
{
flagf=1;
printf(" PowerLaw dN/dr = r** -a \n Exponent = ");
scanf("%s",s);
a1=atof(s);
};
return pow(x,-a1);}
//----------------------------------------
double McCl(double x)
{
char s[20];
if (!flagf)
{
flagf=1;
printf(" McClatchey dN/dr = const or r** -a \n Rmc=");
scanf("%s",s);
a1=atof(s);
printf(" Pokazatel' =");
scanf("%s",s);
a2=atof(s);
};
if (x>a1)
return pow(a1/x,a2);
return 1;
}
//----------------------------------------
double ModGamma(double x)
{
char s[20];
double y;

if (!flagf)
{
flagf=1;
printf(" dN/dr = r**a*exp(-b*r**g) \n a=");
scanf("%s",s);
a1=atof(s);
printf(" b =");
scanf("%s",s);
a2=atof(s);
printf(" g =");
scanf("%s",s);
a3=atof(s);
};

y= log(x);
y= exp(a1*x)*exp(-a2*exp(a3*x));
return y;
}
//----------------------------------------
double LogN(double x)
{
char s[20];
double y;
if (!flagf)
{
flagf=1;
printf(" dN/d_ln(r)=exp(-0.5*((ln(r)-ln(rm))/ln(sigma))^2) \n rm=");
scanf("%s",s);
a1=atof(s);
a1=log(a1);
printf(" sigma =");
scanf("%s",s);
a2=atof(s);
a2=log(a2);
};
y= (log(x)-a1)/a2;
if (fabs(y)>41)
return 0;
return exp( - 0.5*y*y)/x;
}
//----------------------------------------
double Gauss(double x)
{
char s[20];

double y;
if (!flagf)
{
flagf=1;
printf(" dN/dr = exp(-0.5*((r-rm)/sigma)^2) \n rm=");
scanf("%s",s);
a1=atof(s);
printf(" sigma =");
scanf("%s",s);
a2=atof(s);
};
y= (x-a1)/a2;
if (fabs(y)>41)
return 0;
return exp( - 0.5*y*y)/x;
}
//--------------------------------------------------------------
void Legendre(int nn, double x,double *Pl,double *Pm,double *Pn)
{ int i;
double ni;
*Pl = 1;
*Pm = 1;
*Pn = x;
for (i=2;i<=nn;i++)
{ ni=1.0/i;
*Pl=*Pm;
*Pm=*Pn;
*Pn= (2-ni) * x * (*Pm) - (1-ni) * (*Pl); } }
//-----------------------------------------------------------
double FirstGues(int n,int z)
{ return cos(Pi*(z-0.25)/(n+0.5)); }
//---------------------------------}
double Newton(int n,double y)
{
double Pl,Pm,Pn;
Legendre(n,y,&Pl,&Pm,&Pn);
return y-(1-y*y)*Pn/(n*(Pm-y*Pn));
}
//--------------------------------------------------------------
void LobHalf( int n, HalfAngleType muu,HalfAngleType vess)
{
int i,k,m;
double Pl,Pm,Pn,lxd,yd;
m=n+n;
for (i=1;i<=n;i++)
{ lxd= 12686478.0;
yd=cos(Pi*(i-.25)/(m+.5));
k=0;
while ((fabs(yd-lxd)>1e-15)&&(k++< 0x0fff ))
{ lxd = yd;
yd = Newton(m,yd); };
muu[i]=yd;
Legendre(m,yd,&Pl,&Pm,&Pn);
lxd=Pm*m;
vess[i]= 2*(1-yd*yd)/(lxd*lxd);
};
muu[0]=1.0;
vess[0]= 0.0;
}
//---------------------------------------------
long int MaxNumb(double x)
{
double y;
y= exp(log(x)/3.0);
y= x + kl*y + 2;
return (((long int)(y+1)) >> 1) << 1;}

//----------------------------------------------
void GetInd(int nh, double mno,
HalfAngleType r1e,HalfAngleType i1e,
HalfAngleType r2e,HalfAngleType i2e,
HalfAngleType r1o,HalfAngleType i1o,
HalfAngleType r2o,HalfAngleType i2o,
HalfAngleType xf, HalfAngleType xb)
{
int i;
double r1,r2,i1,i2;
for (i=0;i<nh;i++)
{ r1=r1e[i]+r1o[i];
i1=i1e[i]+i1o[i];
r2=r2e[i]+r2o[i];
i2=i2e[i]+i2o[i];
xf[i]= mno*(r1*r1+i1*i1+r2*r2+i2*i2);
r1=r1e[i]-r1o[i];
i1=i1e[i]-i1o[i];
r2=r2e[i]-r2o[i];
i2=i2e[i]-i2o[i];
xb[i]= mno*(r1*r1+i1*i1+r2*r2+i2*i2); }
;}
//--------------------------------------
void ReperPsi(long int ns, int sz,int d_n,double x,ReperType Psi)
{
int i,j,nfr,nost;
long int n;
long double p2,p1,p0,xx;
nfr=ns/sz;
nost= ns % sz;
n= ((ns+d_n) << 1) +5;
xx=1/x;
p2=1.0;
p1=0.0;
for (i=1;i<=d_n;i++)
{ n-=2;
p0=n*xx*p1-p2;
p2=p1;
p1=p0;
}
if (nost)
{ Psi[nfr][1]=p2;
Psi[nfr][0]=p1;
for (i=0;i<nost;i++)
{ n-=2;
p0=n*xx*p1-p2;
p2=p1;
p1=p0;
};
} ;
j=nfr;
while (j--!=0)
{ Psi[j][1]=p2;
Psi[j][0]=p1;
for (i=0;i<sz;i++)
{ n-=2;
p0=n*xx*p1-p2;
p2=p1;
p1=p0;
};
};
p0= 3.0*xx*p1 - p2;
p2= p1;
p1= p0;
if (fabs(p2) > fabs(p1))
p0= (sin(x)/x-cos(x))/p2;
else p0= sin(x)/p1;
if (nost) nfr++;
for (i= 0;i<nfr;i++)
{
Psi[i][0]= p0*Psi[i][0];
Psi[i][1]= p0*Psi[i][1];
};
};
//-------------------------------
void NextPsi(long int n1,int sz,double x,double p1,double p2,RBType8 Psi)
{
long int n;
int i;
long double xx,p0;
i=sz;
n= ((n1+sz)<<1)+3;
xx=1/x;
while (i--!=0)
{
n-=2;
p0=n*xx*p1-p2;
p2= p1;
p1= p0;
Psi[i]= p0;
};
};
//-----------------------------------------------------------------------
void ReperD(long int nmax,long int d_n,int sz,
double rx,double ix,ReperType RDn)
{
long int i,nn;
long double id,rd,xni,xnr,invix,invrx,ssq;
int nfr,nost,j;
ssq= 1/(rx*rx + ix*ix);
invrx= rx * ssq;
invix= - ix * ssq;
nn= d_n + nmax + 2;
rd= nn*invrx;
id= nn*invix;
for (i=0;i<d_n;i++)
{
nn--;
xnr= nn*invrx;
xni= nn*invix;
rd+= xnr;
id+= xni;
ssq= 1/(rd*rd+id*id);
rd= xnr - rd*ssq;
id= xni + id*ssq;
};
nfr= nmax/sz;
nost= nmax%sz;
if (nost)
{
RDn[nfr][0]= rd;
RDn[nfr][1]= id;
if (!nfr) return;
for (i=0;i<nost;i++)
{
nn--;
xnr= nn*invrx;
xni= nn*invix;
rd+= xnr;
id+= xni;
ssq= 1/(rd*rd+id*id);
rd= xnr - rd*ssq;
id= xni + id*ssq;
};
};
for (j=nfr-1;j>0;j--)
{
RDn[j][0]= rd;
RDn[j][1]= id;
for (i=0;i<sz;i++)
{
nn--;
xnr= nn*invrx;
xni= nn*invix;
rd+= xnr;
id+= xni;
ssq= 1/(rd*rd+id*id);
rd= xnr - rd*ssq;
id= xni + id*ssq;
};
};
RDn[0][0]= rd;
RDn[0][1]= id;
};
//-----------------------------------------------------------------------}
void NextD(long int n1,int sz,double rx,double ix,
long double dlr,long double dli,QwType Dn)
{
long int nn;
long double id,rd,xni,xnr,invix,invrx,ssq;
int j;
ssq= 1/(rx*rx + ix*ix);
invrx= rx * ssq;
invix= - ix * ssq;
nn= n1 + sz + 1;
rd= dlr;
id= dli;
for (j=sz-1;j>=0;j--)
{
nn--;
xnr= nn*invrx;
xni= nn*invix;
rd+= xnr;
id+= xni;
ssq= 1/(rd*rd+id*id);
rd= xnr - rd*ssq;
id= xni + id*ssq;
Dn[j][0]= rd;
Dn[j][1]= id;
};
};
//----------------------------------------
void AddEven(int sz,QwType abr,QwType abi,QwType PiTau,
double *i1r,double *i1i,double *i2r,double *i2i)
{
long double x,s1,s2,s3,s4;
int i,j;

s1= *i1r;
s2= *i1i;
s3= *i2r;
s4= *i2i;
j= 0;
for (i= 0; i< (sz >>1); i++)
{
x = PiTau[j][0];
s1+= x * abr[j][0];
s2+= x * abi[j][0];
s3+= x * abr[j][1];
s4+= x * abi[j][1];
j++;
x= PiTau[j][1];
s1+= x * abr[j][1];
s2+= x * abi[j][1];
s3+= x * abr[j][0];
s4+= x * abi[j][0];
j++;
};
if (sz%2)
{
x= PiTau[j][0];
s1+= x * abr[j][0];
s2+= x * abi[j][0];
s3+= x * abr[j][1];
s4+= x * abi[j][1];
};
*i1r= s1;
*i1i= s2;
*i2r= s3;
*i2i= s4;
};
//----------------------------------------
void AddOdd(int sz,QwType abr,QwType abi,QwType PiTau,
double *i1r,double *i1i,double *i2r,double *i2i)
{
long double x,s1,s2,s3,s4;
int i,j;

s1= *i1r;
s2= *i1i;
s3= *i2r;
s4= *i2i;
j= 0;
for (i= 0; i< (sz >>1); i++)
{
x = PiTau[j][1];
s1+= x * abr[j][1];
s2+= x * abi[j][1];
s3+= x * abr[j][0];
s4+= x * abi[j][0];
j++;
x= PiTau[j][0];
s1+= x * abr[j][0];
s2+= x * abi[j][0];
s3+= x * abr[j][1];
s4+= x * abi[j][1];
j++;
};
if (sz%2)
{
x= PiTau[j][1];
s1+= x * abr[j][1];
s2+= x * abi[j][1];
s3+= x * abr[j][0];
s4+= x * abi[j][0];
};
*i1r= s1;
*i1i= s2;
*i2r= s3;
*i2i= s4;
};

//-----------------------------------------------------
void GetPiTau(int sz,long int n,double mu,
double *Pim1, double *Pi0,QwType PT)
{
int i;
long double x,y;
long int nm1,tnm1,np1;
x= mu * (*Pi0);
y= n * (*Pim1);
nm1= n - 1;
if (!nm1) nm1= 1;
np1 = n;
tnm1= (n <<1) - 1;
for (i=0;i<sz;i++)
{
n= np1;
np1++;
*Pim1= *Pi0;
*Pi0= ( tnm1 * x - y )/nm1;
nm1= n;
tnm1= tnm1 + 2;
y= np1 * (*Pim1);
x= mu * (*Pi0);
PT[i][0]= *Pi0;
PT[i][1]= x * n - y;

};
};
//-------------------
void GetM(double br,double bi,double ps0,double ps1,
double x0,double x1,double *cr,double *ci)
{
long double rch,ich,rzn,izn,z;
rch= br * ps1 - ps0;
ich= bi * ps1;
rzn= rch + bi * x1;
izn= ich - br * x1 + x0;
z= rzn*rzn + izn*izn;
*cr= ( rch * rzn + ich * izn) /z;
*ci= ( rzn * ich - rch * izn) /z;
};
//-----------------------------------------------------------------
void BigPR(double m1,double m2,double x,int nh,
double *QE,double *QS,double *QG,double *QB,
HalfAngleType i1re,HalfAngleType i1ie,
HalfAngleType i2re,HalfAngleType i2ie,
HalfAngleType i1ro,HalfAngleType i1io,
HalfAngleType i2ro,HalfAngleType i2io)
{
double m1x,m2x,xx;
long int n_dd,L0,L,n,nmax;
int nframe,nost,sz,n_dp;
long double z,y,cf;
double sni,snr,Qen,d1i,d1r,
ani,anr,bni,bnr,
nix,inmi,inmr,mi,mr,p1,p0,ix,xi1,xi0,xim;
int i,j,k;
ReperType rdn;
ReperType rpsi;
RBType8 psi;
QwType D;
QwType PiTauC;
QwType abr;
QwType abi;

long int tnp1= 1;
double sqs= 0;
double sqe= 0;
double sbr= 0;
double sbi= 0;
double sgQs= 0;
double am1i= 0;
double am1r= 0;
double bm1i= 0;
double bm1r= 0;

double pok=1/3.0;

m1x= m1 * x;
m2x= m2 * x;
mr= m1;
mi= m2;
ix= m1*m1+m2*m2;
inmr= m1/ix;
inmi= - m2/ix;
ix= 1/x;
y= pow(x,pok);
nmax= MaxNumb(x);
n_dp= (kl*y)+15;
if (!(int)m1)
n_dd= n_dp;
else
{
y= kl* ( pow(m1*x,pok) - pow(x,pok) );
L0= 2 + ((m1-1)*x + y)+n_dp;
if (m2 > 1e-15)
y= 30 + 20*m2 + m1/m2*(20.3*m1 - 5.8 - 2.4/(m1-0.8));
else y= 1e9;
if (y > 1e9)
L= 1e9;
else L= y;
if (L < L0) n_dd= L;
else n_dd= L0;

};
// printf("y=%le L=%ld L0=%ld N_DD =%ld\n",(double)y,L,L0,n_dd);


sz= ABSize;
nframe= nmax/sz;
if (nframe > NReper) exit(1);
nost= nmax % sz;
for (i=0;i<nh;i++)
{
PiM1[i]= -1;
Pi0[i]= 0;
i1re[i]= 0;
i1ie[i]= 0;
i2re[i]= 0;
i2ie[i]= 0;
i1ro[i]= 0;
i1io[i]= 0;
i2ro[i]= 0;
i2io[i]= 0;
}
xi0= cos(x);
p0= sin(x);
xim= -p0;

ReperPsi(nmax,sz,n_dp,x,rpsi);
ReperD(nmax,n_dd,sz,m1x,m2x,rdn);
L= 1;
// long int np1= 1;
for (i=0;i<nframe;i++)
{
// printf(" i=%d\n",i);
NextPsi(L,sz,x,rpsi[i][0],rpsi[i][1],psi);
NextD(L,sz,m1x,m2x,rdn[i][0],rdn[i][1],D);
L0= L;
for (j= 0;j<sz;j++)
{
xi1= tnp1*ix*xi0 - xim;
xim= xi0;
tnp1+= 2;
cf= (tnp1*1.0/L)/(L+1);
p1= psi[j];
d1r= D[j][0];
d1i= D[j][1];
nix= ix*L;
snr= nix + d1r*inmr - d1i*inmi;
sni= d1r*inmi + d1i*inmr;
GetM(snr,sni,p0,p1,xi0,xi1,&anr,&ani);
snr= nix + d1r*mr - d1i*mi;
sni= d1r*mi + d1i*mr;
GetM(snr,sni,p0,p1,xi0,xi1,&bnr,&bni);
Qen= (anr+bnr)*tnp1;
sqe+= Qen;
sqs+= (anr*anr + ani*ani + bnr*bnr + bni*bni )*tnp1;
sbr = (anr - bnr)*tnp1 - sbr;
sbi = (ani - bni)*tnp1 - sbi;
sgQs+= (L - 1.0/L) * ( anr * am1r + ani * am1i +
bnr * bm1r + bni * bm1i ) +
cf* ( anr * bnr + ani * bni );
Qen/= sqs;
abr[j][0]= anr * cf;
abr[j][1]= bnr * cf;
abi[j][0]= ani * cf;
abi[j][1]= bni * cf;
am1i= ani;
am1r= anr;
bm1i= bni;
bm1r= bnr;
p0= p1;
xi0= xi1;
L++; // ) {; j:= sz; i:= nframe; - loop exit }
};
for (j= 0; j<nh;j++)
{
GetPiTau(sz,L0,mu[j],&PiM1[j],&Pi0[j],PiTauC);
AddEven(sz,abr,abi,PiTauC,
&i1re[j],&i1ie[j],&i2re[j],&i2ie[j]);
AddOdd(sz,abr,abi,PiTauC,
&i1ro[j],&i1io[j],&i2ro[j],&i2io[j]);
};

};
if (nost)
{
i= nframe;
sz= nost;
NextPsi(L,sz,x,rpsi[i][0],rpsi[i][1],psi);
NextD(L,sz,m1x,m2x,rdn[i][0],rdn[i][1],D);
L0= L;
for (j= 0;j<sz;j++)
{
xi1= tnp1*ix*xi0 - xim;
xim= xi0;
tnp1+= 2;
cf= (tnp1*1.0/L)/(L+1);
p1= psi[j];
d1r= D[j][0];
d1i= D[j][1];
nix= ix*L;
snr= nix + d1r*inmr - d1i*inmi;
sni= d1r*inmi + d1i*inmr;
GetM(snr,sni,p0,p1,xi0,xi1,&anr,&ani);
snr= nix + d1r*mr - d1i*mi;
sni= d1r*mi + d1i*mr;
GetM(snr,sni,p0,p1,xi0,xi1,&bnr,&bni);
Qen= (anr+bnr)*tnp1;
// { if L>x then writeln(L:4,' ',qen,' ',qen/sqe);}


sqe+= Qen;
sqs+= (anr*anr + ani*ani + bnr*bnr + bni*bni )*tnp1;
sbr= (anr - bnr)*tnp1 - sbr;
sbi= (ani - bni)*tnp1 - sbi;
sgQs+= (L - 1.0/L) * ( anr * am1r + ani * am1i +
bnr * bm1r + bni * bm1i ) +
cf* ( anr * bnr + ani * bni );
// printf("L=%ld Delta=%e Eps=%e\n",L,Qen,Qen/sqe);
Qen= Qen/sqs;
// { write(qen:20);}
abr[j][0]= anr * cf;
abr[j][1]= bnr * cf;
abi[j][0]= ani * cf;
abi[j][1]= bni * cf;

am1i= ani;
am1r= anr;
bm1i= bni;
bm1r= bnr;
p0= p1;
xi0= xi1;
L++; // {; j:= sz; i:= nframe; - loop exit }
};

for (j= 0; j<nh;j++)
{
GetPiTau(sz,L0,mu[j],&PiM1[j],&Pi0[j],PiTauC);
AddEven(sz,abr,abi,PiTauC,
&i1re[j],&i1ie[j],&i2re[j],&i2ie[j]);
AddOdd(sz,abr,abi,PiTauC,
&i1ro[j],&i1io[j],&i2ro[j],&i2io[j]);
};
};

ix*= ix;
*QE = sqe*2*ix;
*QS= sqs*2*ix;
*QG= sgQs*4*ix/(*QS);
*QB= ix*(sbr*sbr + sbi*sbi);
};

//-------------------------
void PlusStocks(int nfi,double mno,
HalfAngleType i1re,HalfAngleType i1ie,
HalfAngleType i2re,HalfAngleType i2ie,
HalfAngleType i1ro,HalfAngleType i1io,
HalfAngleType i2ro,HalfAngleType i2io,
AngleType p1,AngleType p2,
AngleType p3,AngleType p4)

{
int i,j=nfi,nh=(nfi+1)>>1;
double r1,r2,i1,i2;
for (i=0;i<nh;i++)
{
j--;
r1= i1re[i] + i1ro[i];
i1= i1ie[i] + i1io[i];
r2= i2re[i] + i2ro[i];
i2= i2ie[i] + i2io[i];

p1[i]+= mno * (r1*r1 + i1*i1);
p2[i]+= mno * (r2*r2 + i2*i2);
p3[i]+= mno * (r1*r2 + i1*i2);
p4[i]+= mno * (r1*i2 - r2*i1);

if (j!=i)
{
r1= i1re[i] - i1ro[i];
i1= i1ie[i] - i1io[i];
r2= i2re[i] - i2ro[i];
i2= i2ie[i] - i2io[i];

p1[j]+= mno * (r1*r1 + i1*i1);
p2[j]+= mno * (r2*r2 + i2*i2);
p3[j]+= mno * (r1*r2 + i1*i2);
p4[j]+= mno * (r1*i2 - r2*i1);
};
};
};

//--------------------------------------------
void BarIntegral(long int n,int nfi,double x1,double x2,
double la,double mr,double mi,double n0,fun DFun,
double *Qe,double *Qs,double *Qg,double *Qb,
AngleType p1,AngleType p2,AngleType p3,AngleType p4)
{
int j;
long int i;
double qe1,qs1,qg1,qb1,f,x,d_x,r,df,mno;
double SBAR = 0;
double SQE = 0;
double SQS = 0;
double SQG = 0;
double SQB = 0;
int nh;
HalfAngleType i1re;
HalfAngleType i1ie;
HalfAngleType i2re;
HalfAngleType i2ie;
HalfAngleType i1ro;
HalfAngleType i1io;
HalfAngleType i2ro;
HalfAngleType i2io;

for (j=0;j<nfi;j++)
{ p1[j]=0;
p2[j]=0;
p3[j]=0;
p4[j]=0;
};

nh=(nfi+1) >> 1;
f= 0.5*la/Pi;
d_x= (x2-x1)/n;
x= x1 + d_x*0.5;
for (i=0;i<n;i++)
{
if (show)
if (!(i%show))
printf("%7d ",i);

x2=x*x;
r= x*f;
df= DFun(r);
BigPR(mr,mi,x,nh,
&qe1,&qs1,&qg1,&qb1,
i1re,i1ie,i2re,i2ie,
i1ro,i1io,i2ro,i2io);

// mno= 2*df/(qs1*x2);
mno=df*d_x;
SQE+= qe1 * mno*x2;
SQS+= qs1 * mno*x2;
SQG+= qs1 * qg1*mno*x2;
SQB+= qb1 * mno*x2;
SBAR+= mno*f;
PlusStocks(nfi,mno,
i1re, i1ie,i2re, i2ie,
i1ro, i1io,i2ro, i2io,
p1, p2, p3, p4);
x+= d_x;
};
if (show)
printf("\n");

mno= 4/SQS;
for (j=0;j<nfi;j++)
{ p1[j]*=mno;
p2[j]*=mno;
p3[j]*=mno;
p4[j]*=mno;
};
mno= la*la*la/SBAR;
*Qe= mno * SQE;
*Qs= mno * SQS;
*Qg= SQG/SQS;
*Qb= mno * SQB;

x2= 0.125/Pi/Pi*n0*1e-12;
*Qe*= x2;
*Qs*= x2;

};
//--------------------------------------------------------
void KbdImp(int *ifun)
{
char ss[20];
printf(" NFiPoint=");
scanf("%s",ss);
nfi=atoi(ss);
if (nfi>NFiMax ) exit(1);


printf(" wavelength[mkm]=");
scanf("%s",ss);
la=atof(ss);

printf(" r_min[mkm]=");
scanf("%s",ss);
r1=atof(ss);

printf(" r_max[mkm]=");
scanf("%s",ss);
r2=atof(ss);


printf(" nrefr_RE=");
scanf("%s",ss);
mr=atof(ss);

printf(" nrefr_Im=");
scanf("%s",ss);
mi=atof(ss);

printf(" Consentr[mg/m^3]=");

scanf("%s",ss);
n0=atof(ss);

printf(" N_IntPoint=");
scanf("%s",ss);
np=atoi(ss);

printf(" 0: disfun=Yunge\n1: disfun=McCl\n2: disfun=ModGamma\n3: disfun=LogN\n4: disfun=Gauss\ndefault : exit(1)\n");

printf(" N_Fun_Distr [0-4]=");
scanf("%s",ss);
*ifun=atoi(ss);
};
//----------------------------------------------------------
void FInp(int *ifun)
{
char ss[90];
double mm1;
int i;
flagf=1;
fgets(ss,sizeof(ss),otx);
fscanf(otx,"%d",&nfi);
fgets(ss,sizeof(ss),otx);

mm1=1;
fscanf(otx,"%lf",&mr);
fscanf(otx,"%lf",&mi);
fgets(ss,sizeof(ss),otx);
if ((fabs(mr)+fabs(mi))!=0)
for (i=0;i<2;i++)
fgets(ss,sizeof(ss),otx);
else
{
fscanf(otx,"%lf",&mm1);
fgets(ss,sizeof(ss),otx);
fscanf(otx,"%lf",&mr);
fscanf(otx,"%lf",&mi);
fgets(ss,sizeof(ss),otx);
mr=mr/mm1;
mi=mi/mm1;
};
fscanf(otx,"%lf",&la);
fgets(ss,sizeof(ss),otx);
la=la/mm1;
fscanf(otx,"%d",ifun);
fgets(ss,sizeof(ss),otx);

fscanf(otx,"%lf",&a1);
fgets(ss,sizeof(ss),otx);
fscanf(otx,"%lf",&a2);
fgets(ss,sizeof(ss),otx);
fscanf(otx,"%lf",&a3);
fgets(ss,sizeof(ss),otx);
fscanf(otx,"%lf",&r1);
fgets(ss,sizeof(ss),otx);
fscanf(otx,"%lf",&r2);
fgets(ss,sizeof(ss),otx);
fscanf(otx,"%ld",&np);
fgets(ss,sizeof(ss),otx);
fscanf(otx,"%lf",&n0);
fgets(ss,sizeof(ss),otx);
fscanf(otx,"%d",&show);

};
//-----------------------------------------------------------
int main(void)

{
double x1,x2,qe,qs,qg,qb;
int i;
double (*disfun)(double);
double x;
int nh;
AngleType p1,p2,p3,p4;

printf("\n");


if ( ( otx=fopen("input.txt","r")) != NULL)
{
FInp(&i);
fclose(otx);
}
else KbdImp(&i);

switch (i)
{ case 0: disfun=Yunge;
sdf="Yunge\0";
break;
case 1: disfun=McCl;
sdf="McCl\0";
break;
case 2: disfun=ModGamma;
sdf="ModGamma\0";
break;
case 3: disfun=LogN;
sdf="LogN\0";
break;
case 4: disfun=Gauss;
sdf="Gauss\0";
break;
default : exit(1);
};

nh=(nfi+1)>>1;
if (nh > (HalfFiMax) )
{ printf(" New_Fi=%d\n",NFiMax);
nh=(nfi+1)>>1;
};
for (i=0;i<nh;i++)
{
x=i*Pi/(nfi-1);
mu[i]=cos(x);
};


x1=2*Pi*r1/la;
x2=2*Pi*r2/la;


BarIntegral(np,nfi,x1,x2,
la,mr,mi,n0,//Yunge, //
disfun,
&qe,&qs,&qg,&qb,
p1,p2,p3,p4);
printf(" Distr Fun = %s\n",sdf);

printf(" QE=%le\n QS=%le\n QG=%le\n QB=%le\n",qe,qs,qg,qb);


otx=fopen("ind_ci.prn","w");
fprintf(otx," N_Integr_P=%ld r1=%lf r2=%lf\n",np,r1,r2);
fprintf(otx," La=%10.7lfmkm m1=%lf m2=%lf n0=%20le\n",la,mr,mi,n0);
fprintf(otx,"Distr Fun = %s %lf %lf %lf\n",sdf,a1,a2,a3);
fprintf(otx," QE=%le\n QS=%le\n QG=%le\n QB=%le\n",qe,qs,qg,qb);
fprintf(otx," Angle I0 S1 S1 S2 S4\n");
for (i=0;i<nfi;i++)
{
x=i*180.0/(nfi-1);
fprintf(otx," %7.2lf %le %le %le %le %le\n",x,(p1[i]+p2[i])/2,
p1[i],p2[i],p3[i],p4[i]);

};


fclose(otx);

// printf("I=%d\n",i);
return 0;
}