Программа расчетов характеристик тумана при рассеянии его искусственными каплями воды
//...................................................................................................................................................
#include <iostream>
#include <stdio.h>
#include <conio.h>
#include <math.h>
#include <stdlib.h>
#include <process.h>
#include <iostream>
#include <fstream>

FILE *AZ;
FILE *REZ;
using namespace std;
double fun(double x, double y);
double my_read(FILE *AA);
int my_error;

double r1[3];
double m1[3];
double dm1[3];
double n10[3];
double f1[3];
double nn1[3];
double V1=0;
double delr1;
double n1,fn1;
double n2;
double r2;
double m2;
double delm;
double V2;
double E[3];
double f2;
double w1,w2;
double l,l0;
double beta1;
double H,Q1,Q2,Q3,Q4,y;
double t, xxx;

//...................................................................................................................................................
int main()
{
setlocale(0,"Russian");
if (( AZ=fopen("az.dat","r")) == NULL){
 printf("error open file AZ.");
 exit(1);
 };
n2=(double)my_read(AZ);
r2=(double)my_read(AZ);
V2=(double)my_read(AZ);
E[0]=(double)my_read(AZ);
E[1]=(double)my_read(AZ);
E[2]=(double)my_read(AZ);
r1[0]=(double)my_read(AZ);
r1[1]=(double)my_read(AZ);
r1[2]=(double)my_read(AZ);
n10[0]=(double)my_read(AZ);
n10[1]=(double)my_read(AZ);
n10[2]=(double)my_read(AZ);
fclose(AZ);
if ((REZ=fopen("rez.dat","w")) ==NULL) {
    printf("Error open file rez.");
    exit(1);
    };
fprintf(REZ,"Исходные данные для крупных капель\n");
fprintf(REZ,"n2 - концентрация (м-3)           %.2e \n",n2);
m2=4.*3.14*pow(r2,3)*1000000/3.;
w2=n2*m2;
fprintf(REZ,"W2 - водность  (г/м3)             %.2e \n",w2);
fprintf(REZ,"r2 - радиус   (мм)                %.2e \n", 1000*r2);
fprintf(REZ,"V2 - скорость падения  (м/с)      %.2e \n",V2);
delm=m2-(4.*3.14*pow((r2-0.0001),3)*1000000/3.);
f2=n2/delm;
for (int i=0; i<3; i++)
{
m1[i]=4.*3.14*pow(r1[i],3)*1000000/3.;
}
delr1=r1[0]+(r1[1]-r1[0]);  //r1[1]-r1[0];
dm1[0]=m1[0]+4.*3.14*pow(100*delr1,3)/3.;
dm1[1]=m1[1]-m1[0];
dm1[2]=m1[2]-m1[1];
n1=0;
for (int i=0; i<3; i++)
{
f1[i]=n10[i]/dm1[i];
n1=n1+n10[i];
}
t=0;
fn1=0;
w1=0;
for (int ii=0; ii<3; ii++)
{
fn1=fn1+n10[ii]*pow(r1[ii],3);
w1=w1+n10[ii]*m1[ii];
}
l0=0.62/fn1;
fprintf(REZ,"Исходные данные для мелких капель\n");
fprintf(REZ,"n1 - концентрация    (м3)        %.2e \n",n1);
fprintf(REZ,"W1 - водность        (г/м3)      %.2e \n",w1);
fprintf(REZ,"r1 - радиус   (мкм)              %.2e, %.2e, %.2e \n",1.e+06*r1[0],1.e+06*r1[1],1.e+06*r1[2]);
fprintf(REZ,"\n Характеристики тумана для различных значений концентрации тумана и осадков\n");
fprintf(REZ," Время         Концентрация     Видимость        Водность\n");
fprintf(REZ," t,сек         n1,м3            МДВ,м            W1,г/м3\n");
fprintf(REZ,"%.2e       %.2e       %.2e          %.2e \n",t,n1,l0,w1);
H=1;
for(int k=5; k<121; k+=5)
{
w1=0;
fn1=0;
n1=0;
for(int ii=0; ii<3; ii++)
{
//...................................................................................................................................................
y=f1[ii];
//...................................................................................................................................................
t=0;
for(int i=0; i<k; i++)
{
Q1=fun(t,y);
Q2=fun(t+H/2.,y+H*Q1/2.);
Q3=fun(t+H/2.,y+H*Q2/2.);
Q4=fun(t+H,y+H*Q3);
y=y+(H/6.)*(Q1-2*Q2+2*Q3+Q4);
t=t+H;
}
nn1[ii]=y*(dm1[ii]);
w1=w1+nn1[ii]*m1[ii];
fn1=fn1+nn1[ii]*pow(r1[ii],2);
n1=n1+nn1[ii];
l=0.62/fn1;
fprintf(REZ,"%.2e       %.2e       %.2e          %.2e \n",t,n1,l,w1);
}
fprintf(REZ,"\n");
fclose (REZ);
system("Pause");
return 0;
}
double fun(double x, double y)
{
double f;
f=-y*beta1;
return(f);
}
double my_read(FILE *AA)
{
static char b[20];
int J;
char c=0;
J=0;
while( (c>'9' || c<'0') && c!='+' &&  c!='-' && c!='E' && c!='.' && !feof(AA) )
c=fgetc(AA);
do{
b[J]=c;
J++;
c=fgetc(AA);
}while( (c<='9' && c>='0') || c=='+' ||  c=='-' || c=='E' || c=='.');
b[J]=0;
if(feof(AA)) my_error=1;
else my_error=0;
return atof(b);
}