/**************************************************************/
/*                                                            */
/* KP                                                         */
/*                                                            */
/* Program to calculate warhead kill probabilities against    */
/* hardened targets. Formulas used were taken from:           */
/*      Tsipis, Kosta. 1983. ARSENAL: UNDERSTANDING           */
/*      WEAPONS IN THE NUCLEAR AGE. New York. Simon and       */
/*      Schuster. Appendix H. Derivation of Kill Probability  */
/*      Equation (pp. 305-308).                               */
/* Anyone may copy or modify this program, as long as the     */
/* result is not used for any profit making purpose.          */
/* Richard J. Stoll, Summer 1987.                             */
/*                                                            */
/**************************************************************/
#include <stdio.h>
#include "math.h"
#include "ctype.h"
#include "conio.h"
#define TRUE 1
#define FALSE 0
main()
{
double Pkp;                        /* kill prob 1 warhead, rel = 1 */
double Pk;                         /* kill prob for 1 warhead, rel < 1 */
double Pk2s,Pk2d;                  /* kp for 2; same, different launchers */
double lr;                         /* lethal radius */
double yield;                      /* yield of warhead */
double EMT;                        /* Equivalent megatonnage of warhead */
double CEP;                        /* CEP of weapon */
double rel;                        /* reliability of weapon */
double H;                          /* hardness of target */
double kv(double y, double cep);   /* Tsipis K value */
double K;                          /* K value */
double fh(double h);               /* hardness function (for lethal radius) */
int more;                          /* T-F: more calcs to do */
FILE *prter;                       /* to use printer */
int Yprt;                          /* T-F: use printer */
int ans;                           /* answer for prompt */


more = TRUE;     /* initialize for more data */
printf("\nDo you want results output to printer (Y/N)? ");
ans = getche();
switch(toupper(ans)){
     case 'Y' : /* get ready for printer */
              prter = fopen("PRN","wt");
              Yprt = TRUE;
              printf("\n**** MAKE SURE PRINTER IS TURNED ON ****");
              printf("\nPress Return when ready");
              getche();		
              break;
     case 'N' : /* no printer */
              Yprt = FALSE;
              break;
}

while(more){  /* more calculations */
    printf("\n\nPARAMETERS");
    printf("\n\tYield of warhead (MT)    = "); scanf("%lf",&yield);
    printf("\tCEP of weapon (nm)       = "); scanf("%lf",&CEP);
    printf("\tReliability of weapon    = "); scanf("%lf",&rel);
    printf("\tHardness of target (lbs) = "); scanf("%lf",&H);

    if(Yprt){
          fprintf(prter,"\n\nPARAMETERS");
          fprintf(prter,"\n\tYield of warhead (MT)    = %9.3lf",yield);
          fprintf(prter,"\n\tCEP of weapon (nm)       = %9.3lf",CEP);
          fprintf(prter,"\n\tReliability of weapon    = %9.3lf",rel);
          fprintf(prter,"\n\tHardness of target (lbs) = %9.3lf",H);
    }

    /* LETHAL RADIUS */
    if(H > 1000)
        lr = pow(yield,0.33)/(pow(H,0.33)*pow(fh(H),0.33));
    else  /* don't use fh() if H too small */
        lr = pow((16.*(yield/H)),0.33);

    /* K VALUE */
    K = kv(yield,CEP);


    /* EMT of warhead */
    EMT = pow(yield,(2.0/3.0));

    /* Pk VALUE, rel = 1 */
    Pkp = 1.0 - exp(-K/(0.22*pow(H,0.67)));

    /* Pk VALUE, rel < 1 */
    Pk = rel*Pkp;

    /* Pk2, one launcher */
    Pk2s = rel*(1.0 - (1.0 - Pkp)*(1.0 - Pkp));

    /* Pk2, two launchers */
    Pk2d = 1.0 - (1.0 - rel*Pkp)*(1.0 - rel*Pkp);

    printf("\n\nCALCULATIONS");
    printf("\n\tLethal Radius     = %8.3lf",lr);
    printf("\n\tEMT of warhead    = %8.3lf",EMT);
    printf("\n\tValue of K        = %8.3lf",K);
    printf("\n\tPk,  rel of %4.2lf  = %8.3lf",1.0,Pkp);
    if(rel < 1.0)
       printf("\n\tPk,  rel of %4.2lf  = %8.3lf",rel,Pk);
    printf("\n\tPk2, 1 launcher   = %8.3lf",Pk2s);
    printf("\n\tPk2, 2 launchers  = %8.3lf",Pk2d);

    if(Yprt){
         fprintf(prter,"\n\nCALCULATIONS");
	 fprintf(prter,"\n\tLethal Radius     = %8.3lf",lr);
	 fprintf(prter,"\n\tEMT of warhead    = %8.3lf",EMT);
         fprintf(prter,"\n\tValue of K        = %8.3lf",K);
         fprintf(prter,"\n\tPk,  rel of %4.2lf  = %8.3lf",1.0,Pkp);
         if(rel < 1.0)
            fprintf(prter,"\n\tPk,  rel of %4.2lf  = %8.3lf",rel,Pk);
         fprintf(prter,"\n\tPk2, 1 launcher   = %8.3lf",Pk2s);
         fprintf(prter,"\n\tPk2, 2 launchers  = %8.3lf",Pk2d);
         fprintf(prter,"\n");
    }

    printf("\n\nMore calculations (Y/N)? ");
    ans = getche();
    switch(toupper(ans)){
       case 'Y' :
             more = TRUE;
             break;
       case 'N' :
             more = FALSE;
             break;
    }

} /* more calculations */

} /* end of main() */

double kv(double y, double cep)
/**************************************************************/
/*                                                            */
/* kv                                                         */
/*                                                            */
/* Calculates Tsipis K value, a relative measure of the       */
/* lethality of a warhead against a hardened target. It is    */
/* relative in the sense that the K values for 2 warheads can */
/* be compared, and the larger K is the better warhead, but   */
/* they may both be lousy or both be great against a target   */
/* with a particular K value.                                 */
/*                                                            */
/**************************************************************/
{
double k;    /* holds k value */

k = pow(y,0.67)/(cep*cep);

return(k);

} /* *** end of kv() *** */

double fh(double h)
/**************************************************************/
/*                                                            */
/* fh                                                         */
/*                                                            */
/* Strange little function that appears in the calculation of */
/* the lethal radius. See Tsipis, page 306.                   */
/*                                                            */
/**************************************************************/
{
double tmp;     /* holds function value */

tmp = -0.19*(1.0/h) - 0.23*(pow(h,-0.5)) + 0.068;

return(tmp);

} /* *** end of fh() *** */