/*
 * @(#) smt1_4.c - Program for calculation of the hitting factor of
 * an air bomb (Task 1_4 at the Special Military Training course).
 * (c) 1996 Ivan Maidanski <ivmai@chat.ru> http://ivmai.chat.ru
 * Freeware program source. All rights reserved.
 **
 * Language: ANSI C
 * Tested with: Borland C++ v3.1
 * Last modified: 1996-12-12 19:45:00 GMT+03:00
 */

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

#define K 12

int rnd_func(int x)
{
 return (x*5%(1<<K));
}

double f_func(double x, double y, double z,
              double mx, double my, double mz,
              double sigmax, double sigmay, double sigmaz)
{
 x=(x-mx)/sigmax;
 y=(y-my)/sigmay;
 z=(z-mz)/sigmaz;
 return (exp(-(x*x+y*y+z*z)/2.0)*M_1_PI*M_1_SQRTPI/
         (sigmax*sigmay*sigmaz*sqrt(8.0)));
}

double g_func(double x, double y, double z,
              double a, double b, double k,
              double r_phi, double r_1,
              double tan_phi_1, double tan_phi_2)
{
 x=x*x+y*y;
 y=sqrt(x);
 a=(a+r_phi)*(a+r_phi);
 if (x<a && fabs(z)<b+r_phi)
  return 1.0;
 if (y<r_1 && y*tan_phi_1<z+b && y*tan_phi_2>z-b)
  return (exp((a-x)/k));
 return 0.0;
}

int main()
{
 int i,n, seed;
 double x,y,z,mz,sigmax,sigmay,sigmaz,l;
 double r_phi,r_1,phi_1,phi_2,a,b,k;
 double c_1,c_2;
 printf("Input SigmaX (SigmaX>0):\n");
 scanf("%lf",&sigmax);
 printf("Input SigmaY (SigmaY>0):\n");
 scanf("%lf",&sigmay);
 printf("Input C1 (C1<0):\n");
 scanf("%lf",&c_1);
 printf("Input C2 (0<C2<1):\n");
 scanf("%lf",&c_2);
 printf("Input r_phi (r_phi>0):\n");
 scanf("%lf",&r_phi);
 printf("Input r_1 (r_1>0):\n");
 scanf("%lf",&r_1);
 printf("Input phi_1 (rad):\n");
 scanf("%lf",&phi_1);
 printf("Input phi_2 (rad):\n");
 scanf("%lf",&phi_2);
 printf("Input a (a>0):\n");
 scanf("%lf",&a);
 printf("Input b (b>0):\n");
 scanf("%lf",&b);
 printf("Input K (K>0):\n");
 scanf("%lf",&k);
 printf("Input N (10<=N<=1000):\n");
 scanf("%d",&n);
 printf("Input Seed [0..%d]:\n",(1<<(K-3))*3);
 scanf("%d",&seed);
 seed=(seed<<1)+(1<<(K-2))-1;
 phi_1=tan(M_PI_2-phi_1);
 phi_2=tan(M_PI_2-phi_2);
 l=0.0;
 for (i=0;i<n;i++)
 {
  x=((double)(seed=rnd_func(seed))/(1<<K)*2-1)*sigmax*3;
  y=((double)(seed=rnd_func(seed))/(1<<K)*2-1)*sigmay*3;
  z=sqrt(x*x+y*y);
  mz=z*c_1;
  sigmaz=z*c_2;
  z=((double)(seed=rnd_func(seed))/(1<<K)*2-1)*sigmaz*3+mz;
  l+=f_func(x,y,z,0,0,mz,sigmax,sigmay,sigmaz)*sigmaz*
      g_func(x,y,z,a,b,k,r_phi,r_1,phi_1,phi_2);
 }
 printf("\n W=%lf\n\n",216.0*l/n*sigmax*sigmay);
 return 0;
}