IZP Projekt 2

/*
 * Soubor:  proj2.c
 * Datum:   17.11.2008
 * Autor:   David Sabata, xsabat01@stud.fit.vutbr.cz
 * Projekt: Aproximace funkci a slovni ulohy, projekt ź. 2 pro pýedmŘt IZP
 * Popis:   Program provadi aproximaci funkci LN a TAN se zvolenou odchylkou
 *          a dale zpracovava udaje ze dvou radaru, a to konkretne prepoctem
 *          polohy na zrychleni a zrychleni na polohu.   
 */
 
 
// prace se vstupem/vystupem
#include <stdio.h>
 
// obecne funkce jazyka C
#include <stdlib.h>
 
// kvuli funkci strcmp
#include <string.h>
 
// matematicke funkce
#include <math.h>
 
// Vysledky rozboru parametru
enum actions
{
  AHELP,    // Napoveda
  ATAN,     // Tangens
  ALN,      // Logaritmus
  ARADAR1,  // Radar1
  ARADAR2   // Radar2
};
 
 
// Chybove hlasky
const char *EMSG[] =
{
  /* EOK */ 
  "Vse v poradku.\n",
  /* ECLWRONG */
  "Chybne parametry prikazoveho radku!\n",
  /* ENEGEPS */
  "Chybny parametr! Epsilon musi byt kladne cislo!\n",
};
 
// Pro chybove hlasky
enum codes
{
  EOK = 0,  // Zadna chyba
  ECLWRONG, // Chybne parametry
  ENEGEPS,  // Zaporne EPS
};
 
// Nase konstanty
const double IZP_E = 2.7182818284590452354;        // e
const double IZP_PI = 3.14159265358979323846;      // pi
const double IZP_2PI = 6.28318530717958647692;     // 2*pi
const double IZP_PI_2 = 1.57079632679489661923;    // pi/2
const double IZP_PI_4 = 0.78539816339744830962;    // pi/4
 
const double RADAR_TO = 0.5;  // Doba mezi prichozimi signaly radaru
 
#define RADAR_DIMS 3  // Radar pracuje v 3D prostoru: 1bod=3 souradnice 
 
// Struktura pro kazdy jednotlivy smer bodu (draha, rychlost)
typedef struct {
   double s;   // Poloha
   double v;   // Rychlost
} TRadar_dir;
 
// Struktura pro kazdy bod, obsahuje potrebny pocet smerovych struktur 
// a ukazatel na prave pocitany smer, ktery se pouziva pro predavani funkcim
typedef struct {
   TRadar_dir direction[RADAR_DIMS];   // Pole smerovych struktur
   TRadar_dir *p_dir;   // Ukazatel na aktualni smerovou strukturu
   int actDir;    // Index pole aktualni smerove struktury
} TRadar_dot;
 
////////////////////////////////////////////////////////////////////////////////
void printError(int error);
void printHelp();
int loadParams(int argc, char *argv[]);
void loadInput(int action, double eps, TRadar_dot *radar);
double mySin(double x, double eps);
double myCos(double x, double eps);
double myTan(double x, double eps);
double myLn(double x, double eps);
double convLn(double x, double eps);
double radar1(TRadar_dir *rad, double a);
double radar2(TRadar_dir *rad, double s);
////////////////////////////////////////////////////////////////////////////////
/**
 * Hlavni program
 */
int main(int argc, char *argv[])
{
   // Zjisteni provadene akce na zakladu rozboru parametru
   int action = loadParams(argc, argv);
 
 
   double eps;    // presnost nactena z parametru
 
 
   // Deklarace a inicializace radarovych struktur
   TRadar_dot radar;
   if (action == ARADAR1 || action == ARADAR2)
   {      
      // Vynulovani pocatecnich stavu a nastaveni ukazatele na prvni souradnici 
      for (int i = RADAR_DIMS-1; i>=0; i--) {
         radar.actDir = i;
         radar.p_dir = &radar.direction[radar.actDir];
         radar.p_dir->s = 0;
         radar.p_dir->v = 0;      
      }
   }   
 
 
 
   // Pokud nevypisujeme napovedu, budeme nacitat ze vstupu
   if (action!=AHELP)
   {
      // Nezavisle na provadene akci nacteme epsilon, pokud existuje
      // a kontrolujeme jestli je kladne (mohlo by byt i nula, ale toto je
      // taky chyba, protoze by byl vypocet nekonecny)
      if (argc > 2) 
      {
         sscanf(argv[2], "%lf", &eps);
         if (eps <= 0.0) 
            printError(ENEGEPS);
      } 
 
      // Spustime hlavni nacitaci smycku
      loadInput(action, eps, &radar);
 
   } 
   else 
   {  // Vypis napovedy
      printHelp();
   }
 
 
 
   return EXIT_SUCCESS;
}
 
 
 
/**
 * Vytiskne napovedu 
 */
void printHelp()
{   
   printf(
      "Program Aproximace funkci a radary\n"
      "Autor: David Sabata (c)2008\n"
      "Program aproximuje funkce TAN a LN pomoci Taylorovy rady s presnosti \n"
      "zadanou v parametru. Po spusteni nacita ze vstupu libovolne mnozstvi \n"
      "vstupnich hodnot oddelenych jednim nebo vice bilymi znaky. \n"
      "Druha cast programu pocita polohu nebo zrychleni bodu, sledovaneho \n"
      "radarem ve 3D prostoru. Vstupni hodnoty polohy, oddelene libovolnym \n"
      "poctem bilych znaku, prevadi na hodnoty zrychleni a naopak hodnoty \n"
      "zrychleni prevadi na hodnoty polohy.\n"
      "Parametry:\n"
      "\t-h\t\tzobrazeni teto napovedy\n"
      "\t-ln EPS \tprogram bude aproximovat funkci LN s presnosti EPS\n"
      "\t-tan EPS\tprogram bude aproximovat funkci TAN s presnosti EPS\n"
      "\t-radar1 \tprogram bude prepocitavat zrychleni -> poloha\n"
      "\t-radar2 \tprogram bude prepocitavat poloha -> zrychleni\n\n"
      "Hodnoty parametru i vstupni hodnoty programu muzete vkladat jako \n"
      "cela nebo desetinna cisla v beznem nebo exponencialnim tvaru. \n"
      "Vstupni hodnoty pro funkci TAN zadavejte v radianech.\n"      
   );
}
 
 
/**
 * Vytiskne chybu a ukonci program
 * @param code kod chyby nebo stavu programu
 */
void printError(int error)
{   
   fprintf(stderr, "Chyba: %s\n", EMSG[error]);
   exit(EXIT_FAILURE);
}
 
 
/**
 * Nacte parametry a vyhodnoti akci, pripadne chybu
 * @param argc int pocet parametru
 * @param argv char pole parametru
 * @return int akce ktera se ma provest   
 */ 
int loadParams(int argc, char *argv[])
{
   // Napoveda
   if (argc == 2 && strcmp("-h", argv[1]) == 0) 
      return AHELP;
 
   // Tangens
   if (argc == 3 && strcmp("-tan", argv[1]) == 0) 
      return ATAN;
 
   // Logaritmus
   if (argc == 3 && strcmp("-ln", argv[1]) == 0) 
      return ALN;
 
   // Radar 1
   if (argc == 2 && strcmp("-radar1", argv[1]) == 0) 
      return ARADAR1;
 
   // Radar 2
   if (argc == 2 && strcmp("-radar2", argv[1]) == 0) 
      return ARADAR2;
 
   // Cokoliv jineho je chyba, vytiskneme hlasku o nespravnych parametrech,
   // coz zaroven ukonci program. Hodnota v returnu je pouze formalita. 
   printError(ECLWRONG);
 
   return -1;
}
 
 
/**
 * Hlavni smycka pro nacitani ze vstupu a volani pozadovanych funkci podle
 * predaneho parametru
 * @param action int akce definovana parametry
 * @param eps double presnost definovana parametry
 * @param radar TRadar_bod radarova struktura
 * @return void
 */   
void loadInput(int action, double eps, TRadar_dot *radar)
{
   double inp,       // hodnota nactena ze vstupu 
          tmp;       // pomocna promenna pro vypisovani aproximovanych hodnot
 
   int check;        // kontrola - pocet nactenych hodnot fci scanf
 
   // Nacitame dokud je neco na vstupu
   while ((check = scanf("%lf", &inp)) != EOF) 
   {
      // S neplatnym znakem pocitame jako s NAN a preskocime jeho nacitani
      if (check == 0)
      {
         inp = NAN;
         scanf("%*s");            
      }
 
      switch (action) 
      {                
         // Tangens
         case ATAN:               
            tmp = myTan(inp, eps);            
 
            printf("%.10le\n", tmp);
            break;
 
         // Logaritmus
         case ALN:
            // Fce convLn obaluje fci myLn a pridava prevedeni do konvergence
            tmp = convLn(inp, eps);
 
            printf("%.10le\n", tmp);               
            break;
 
         // Radar1
         case ARADAR1:
            // Vypocet a vypis nove pozice
            tmp = radar1(radar->p_dir, inp);
            printf("%.10le\n", tmp);      
 
            // Posunuti citace na dalsi souradnici
            if (radar->actDir == RADAR_DIMS-1)
               radar->actDir = 0;
            else
               radar->actDir++;
 
            // Prenastaveni ukazatele na dalsi souradnici
            radar->p_dir = &radar->direction[radar->actDir];
 
            break;
 
         // Radar2
         case ARADAR2:
            // Vypocet a vypis noveho zrychleni
            tmp = radar2(radar->p_dir, inp);
            printf("%.10le\n", tmp);               
 
            // Posunuti citace na dalsi souradnici
            if (radar->actDir == RADAR_DIMS-1)
               radar->actDir = 0;
            else
               radar->actDir++;
 
            // Prenastaveni ukazatele na dalsi souradnici
            radar->p_dir = &radar->direction[radar->actDir];
 
            break;
      }
   } 
 
}
 
 
 
/**
 * Aproximacni vypocet fce sinus
 * @param x double uhel pro vypocet v radianech
 * @param eps double odchylka pro urceni presnosti
 * @return double sinus daneho uhlu
 */    
double mySin(double x, double eps)
{
   // Napodobeni matematickych funkci - pro INFINITY nedefinovano
   if (x == INFINITY)
   {
      return NAN;
   }
 
   // Obraceni do minusu - nula neobraci, jednicka obraci
   int trans = 0;
 
   // Dostaneme x na prvni periodu
   x = fmod(x, IZP_2PI);   
 
   // Pujdeme jenom do kladne casti osy
   if (x < 0) 
   {
      x = -x;              // fce sin je licha: sin(-x) = -sin(x)
      trans = 1 - trans;   // zapamatujeme si vynasobeni vysledku -1
   }
 
   // Omezime se na <0; PI> - presuneme x a pamatujeme si otoceni do minusu   
   if (x > IZP_PI) 
   {
      x -= IZP_PI;         // x se zapornym sinem presuneme na prvni pulperiodu
      trans = 1 - trans;   // a zapamatujeme si vynasobeni vysledku -1
   }   
 
   // Dostaneme se do prvniho kvadrantu
   if (x > IZP_PI_2)
   {
      x = IZP_PI - x;
   }
 
   double k = 1;        // citac kroku pro iterace, zvysuje se po 2
   double t = x;        // jeden clen posloupnosti + nastaveni prvni hodnoty
   double sum = t;      // celkovy soucet + ulozeni prvniho clenu
   double xx = x * x;   // x na druhou
 
   // Hlavi smycka
   // sin(x) = x - x^3/3! + x^5/5! - x^7/7! + ...
   while (fabs(t) >= eps) {
      k += 2;
      t = (-t * xx) / ( k * (k - 1));
      sum += t;
   }
 
   // Pokud jsme pri prevadeni do prvniho kvadrantu pouzili lichy pocet 
   // otaceni kolem osy (nasobeni -1), vracime vysledek jako zaporny
   if (trans)
      return -sum;
   else
      return sum;   
}
 
 
/**
 * Aproximacni vypocet fce kosinus
 * @param x double uhel pro vypocet v radianech
 * @param eps double odchylka pro urceni presnosti
 * @return double kosinus daneho uhlu
 */    
double myCos(double x, double eps)
{
   // Napodobeni matematickych funkci - pro INFINITY nedefinovano
   if (x == INFINITY)
   {
      return NAN;
   }
 
   // Obraceni do minusu - nula neobraci, jednicka obraci
   int trans = 0;
 
   // Funkce je suda (cos(x)=cos(-x)), presuneme se do kladne casti osy x
   x = fabs(x);
 
   // Dostaneme x na prvni periodu
   x = fmod(x, IZP_2PI);
 
   // Dostaneme se na prvni pulperiodu (krivka je soumerna podle PI) 
   if (x >= IZP_PI) 
   {
      x = IZP_2PI - x;
   }
 
   // Dostaneme se do prvniho kvadrantu
   if (x > IZP_PI_2) 
   {
      x = -1 * (IZP_PI - x);  // druhy kvadrant je bodove soumerny podle PI/2
      trans = 1 - trans;      // pamatujeme si vynasobeni vysledku -1
   }
 
 
   double k = 0;     // citac kroku pro iterace, zvysuje se po 2
   double t = 1;     // clen posloupnosti + nastaveni prvni hodnoty
   double sum = t;   // celkovy soucet + ulozeni hodnoty prvniho clenu
 
   // Hlavni smycka
   // cos(x) = 1 - x^2/2! + x^4/4! - x^6/6! + ...
   while (fabs(t) >= eps) {
      k += 2;
      t = (-t * x * x) / ( k * (k - 1));
      sum += t;
   }
 
   // Pokud jsme pri prevadeni do prvniho kvadrantu pouzili lichy pocet 
   // otaceni kolem osy (nasobeni -1), vracime vysledek jako zaporny
   if (trans)
      return -sum;
   else
      return sum;
}
 
 
 
/**
 * Funkce pro vypocet tangens; vola mySin a myCos a pdle vysledku vraci bud 
 * jejich podil nebo NAN v pripade nedefinovane hodnoty
 * @param x double hodnota pro kterou pocitame tan (v radianech)
 * @param eps double stanoveni presnosti
 * @return double tangens vstupniho uhlu
 */ 
double myTan(double x, double eps)
{
   double apxSin = mySin(x, eps);
   double apxCos = myCos(x, eps);              
 
   // Overeni definovanosti tan
   if (apxCos == 0)
      return NAN;
   else
      return apxSin/apxCos;   
}
 
 
 
/**
 * Aproximacni vypocet prirozeneho logaritmu (bez prevodu do oblasti 
 * konvergence, o to se stara funkce convLn, ktera po rozkladu vola tuto fci)
 * @param x double cislo pro ktere pocitame logaritmus
 * @param eps double odchylka pro urceni presnosti
 * @return double prirozeny logaritmus cisla
 */    
double myLn(double x, double eps)
{
   // Napodobeni matematickych funkci - pro hodnoty <=0 nedefinovano
   if (x <= 0)
   {
      return NAN;
   }
 
   double t = 1;                 // clen posloupnosti + nastaveni prvniho clenu
   double y = (x - 1) / (x + 1); // y se ve vzorci pouziva jako substituce x
   double yy = y * y;            // y na druhou
   double k = 1;                 // citac pro iterace, zvysuje se po 2
   double sum = t;               // celkovy soucet posloupnosti
 
   // Hlavni smycka
   // ln(x) = 2y*( 1 + y^2/3 + y^4/5 + y^6/7 + ... )
   // kde y = (x-1) / (x+1)
   while (t >= eps) {
      k += 2;
      t = t * yy * (k - 2) / k;
      sum += t;      
   }
 
   return 2 * y * sum;
}
 
 
/**
 * Pocita ln rozlozenim do vice vypoctu v oblasti konvergence
 * @param x double cislo pro ktere pocitame logaritmus
 * @param eps double odchylka pro urceni presnosti
 * @return double prirozeny logaritmus cisla
 */    
double convLn(double x, double eps)
{
   // Kontrola "definovanosti"
   if (x == 0)
      return -INFINITY;
   if (x < 0)
      return NAN;
 
   double ln = 0;
   double convEnd = 2;
   double convStr = 0.25;
 
   // Vsechno vetsi nez horni mez konvergence rozkladame
   if (x > convEnd) 
   {
      // Pomocna hodnota - ln horni meze - tu pak pricitame k vysledku
      double lnConvEnd = myLn(convEnd, eps);
 
      // Dokud se nedostaneme pod horni mez, budeme delit jeji hodnotou
      // a pricitat k vysledku hodnotu jejiho prir. logaritmu
      while (x > convEnd)
      {
         ln += lnConvEnd;
         x = x / convEnd;               
      }    
   }   
 
   // Vsechno mensi nez dolni mez konvergence rozkladame
   if (x < convStr)
   {
      // Pomocna hodnota - ln dolni meze - tu pak odecitame od vysledku
      double lnConvStr = myLn(convStr, eps); 
 
      // Dokud se nedostaneme nad dolni mez, budeme delit jeji hodnotou
      // a pricitat k vysledku hodnotu jejiho prir. logaritmu
      while (x < convStr)
      {
         ln += lnConvStr;
         x = x / convStr;
      }    
   }
 
   // Spocitat zbytek, ktery uz je v oblasti konvergence
   ln += myLn(x, eps);
 
   return ln;
}
 
 
 
/**
 * Prvni radar - prepocitava zrychleni na polohu
 * @param bod TRadar_dot struktura zachycujici stav bodu
 * @param a double nove zrychleni
 * @return double nova poloha   
 */ 
double radar1(TRadar_dir *rad, double a) {
 
   rad->s += 0.5 * RADAR_TO * RADAR_TO * a + rad->v * RADAR_TO;
   rad->v += a * RADAR_TO;      
 
   return rad->s;
}
 
 
 
/**
 * Prvni radar - prepocitava polohu na zrychleni
 * @param bod TRadar_dir struktura zachycujici stav bodu
 * @param s double nova poloha
 * @return double nove zrychleni  
 */ 
double radar2(TRadar_dir *rad, double s) {
 
   double a = 2 * (s - rad->v * RADAR_TO - rad->s) / (RADAR_TO * RADAR_TO);
 
   rad->s = s;
   rad->v += a * RADAR_TO;            
 
   return a;
}