src/Tools/rec3d/ProjectionDirectionMapCreate/src/ProjectionDirectionMapCreate.c

   1 /*
   2 # ProjectionDirectionMapCreate : $Revision$
   3 # $Date$
   4 # Created by $Author$
   5 # Usage : ProjectionDirectionMapCreate
   6 # Attention
   7 #   $Loccker$
   8 #       $State$
   9 #
  10 */
  11 #include <stdio.h>
  12 #include <stdlib.h>
  13 #include <string.h>
  14 #include <math.h>
  15 #define GLOBAL_DECLARATION
  16 #include "../inc/config.h"
  17
  18 #define DEBUG
  19 #include "genUtil.h"
  20 #include "Memory.h"
  21 #include "String.h"
  22 #include "Matrix3D.h"
  23 #define WORDLEN 1024
  24
  25 typedef struct lProjectionDirectionMapCreateInfo {
  26         int LineNum;
  27         char** FileName;
  28         char** RotationalMode;
  29         float* Phi;
  30         float* Theta;
  31         float* Psi;
  32         double* x;
  33         double* y;
  34         double* x2;
  35         double* y2;
  36         int numP;
  37         int numM;
  38         int OutputMode;
  39         int mode;
  40 } lProjectionDirectionMapCreateInfo;
  41
  42 /*
  43 typedef enum lProjectionDirectionMapCreateMode {
  44         a=0,
  45         b=1
  46 } lProjectionDirectionMapCreateMode;
  47 */
  48
  49 void
  50 __InputDataFileListRead(lProjectionDirectionMapCreateInfo* linfo, FILE* fpt, int mode);
  51
  52 void
  53 __RotationalModeChange(lProjectionDirectionMapCreateInfo* linfo, int mode);
  54
  55 void
  56 lProjectionDirectionMapCreate(lProjectionDirectionMapCreateInfo* linfo, int mode);
  57
  58 void
  59 __OutputFileWrite(lProjectionDirectionMapCreateInfo* linf, FILE* fpt, int mode);
  60
  61 int
  62 main(int argc, char* argv[])
  63 {
  64         ProjectionDirectionMapCreateInfo  info;
  65         lProjectionDirectionMapCreateInfo linfo;
  66         int i;
  67
  68         init0(&info);
  69     argCheck(&info, argc, argv);
  70     init1(&info);
  71
  72         DEBUGPRINT("Program Start\n");
  73         if(info.flagOut == 1) linfo.OutputMode = 0;
  74         if(info.flagPs  == 1) linfo.OutputMode = 1;
  75         linfo.mode = info.mode;
  76
  77         fprintf(stderr, "%d\n",linfo.mode);
  78
  79         __InputDataFileListRead(&linfo, info.fptIn, 0);
  80         switch(info.Rmode){
  81                 case 0 :
  82                 break;
  83                 case 1 :
  84                         __RotationalModeChange(&linfo, 0);
  85                 break;
  86         }
  87         lProjectionDirectionMapCreate(&linfo, 0);
  88         switch(linfo.OutputMode){
  89                 case 0 :
  90                         __OutputFileWrite(&linfo, info.fptOut, 0);
  91                         break;
  92                 case 1 :
  93                         fprintf(stderr, "Not supported this mode");
  94                         exit(1);
  95                         break;
  96         }
  97
  98         DEBUGPRINT("Program End\n");
  99         exit(EXIT_SUCCESS);
 100 }
 101
 102 void
 103 additionalUsage()
 104 {
 105         fprintf(stderr, "----- Rotational Mode -----\n");
 106         fprintf(stderr, "Final rotation of projection is only Z .\n");
 107         fprintf(stderr, "\n");
 108         fprintf(stderr, "----- input file format -----\n");
 109         fprintf(stderr, "filename rotationalmode angle angle angle\n");
 110         fprintf(stderr, "   .             .        .     .     .  \n");
 111         fprintf(stderr, "   .             .        .     .     .  \n");
 112         fprintf(stderr, "   .             .        .     .     .  \n");
 113         fprintf(stderr, "   .             .        .     .     .  \n");
 114         fprintf(stderr, "\n");
 115         fprintf(stderr, "----- mode -----\n");
 116         fprintf(stderr, "-m : 0 : Not separate plus and minus of theta. Create one line data\n");
 117         fprintf(stderr, "     1 :     separate plus and minus of theta. Create one line data\n");
 118         fprintf(stderr, "     2 :     separate plus and minus of theta. Create two line data\n");
 119         fprintf(stderr, "     3 : Mollweide expressionn");
 120         fprintf(stderr, "\n");
 121         fprintf(stderr, "----- mode -----\n");
 122         fprintf(stderr, "-Rm : 0 : Not change rotational mode\n");
 123         fprintf(stderr, "      1 :     Change rotational mode to ZONS\n");
 124         fprintf(stderr, "\n");
 125 }
 126
 127 void
 128 __InputDataFileListRead(lProjectionDirectionMapCreateInfo* linfo, FILE* fpt, int mode)
 129 {
 130         int i;
 131         char *s;
 132
 133         fseek(fpt, 0L, SEEK_SET);
 134         s = MemoryAllocate(char, WORDLEN, "in __InputDataFileListRead");
 135         linfo->LineNum = 0;
 136         while(fgets(s, WORDLEN-1, fpt) != NULL){
 137                 linfo->LineNum++;
 138         }
 139         fseek(fpt, 0L, SEEK_SET);
 140
 141         linfo->FileName       = MemoryAllocate(char*, linfo->LineNum, "in __InputDataFileListRead");
 142         linfo->RotationalMode = MemoryAllocate(char*, linfo->LineNum, "in __InputDataFileListRead");
 143         for(i=0; i<linfo->LineNum; i++){
 144                 linfo->FileName[i]       = MemoryAllocate(char, 256, "in __InputDataFileListRead");
 145                 linfo->RotationalMode[i] = MemoryAllocate(char, 5,   "in __InputDataFileListRead");
 146         }
 147         linfo->Phi   = MemoryAllocate(float, linfo->LineNum, "in __InputDataFileListRead");
 148         linfo->Theta = MemoryAllocate(float, linfo->LineNum, "in __InputDataFileListRead");
 149         linfo->Psi   = MemoryAllocate(float, linfo->LineNum, "in __InputDataFileListRead");
 150
 151         for(i=0; i<linfo->LineNum; i++){
 152                 stringGetFromFile(s, "", fpt, stdout, 3);
 153                 stringCopy(linfo->FileName[i],       stringGetNthWord(s, 1, " ,\t"), 255);
 154                 stringCopy(linfo->RotationalMode[i], stringGetNthWord(s, 2, " ,\t"),   5);
 155                 linfo->Phi[i]   = (float)stringGetNthRealData(s, 3, " ,\t")*RADIAN;
 156                 linfo->Theta[i] = (float)stringGetNthRealData(s, 4, " ,\t")*RADIAN;
 157                 linfo->Psi[i]   = (float)stringGetNthRealData(s, 5, " ,\t")*RADIAN;
 158         }
 159         free(s);
 160 }
 161
 162 void
 163 __RotationalModeChange(lProjectionDirectionMapCreateInfo* linfo, int mode)
 164 {
 165         int i, j;
 166         Matrix3D Matrix;
 167         matrix3DParaTypeReal x, y, z;
 168
 169         for(i=0; i<linfo->LineNum; i++){
 170                 matrix3DRotationSetFollowingEulerAngle(Matrix,
 171                                                                                        linfo->RotationalMode[i],
 172                                                                                            linfo->Phi[i],
 173                                                                                            linfo->Theta[i],
 174                                                                                            linfo->Psi[i],
 175                                                                                            MATRIX_3D_MODE_INITIALIZE);
 176                 matrix3DEulerAngleGetFromMatrix3D(Matrix, "ZONS", &x, &y, &z, 0);
 177                 linfo->Phi[i]   = x;
 178                 linfo->Theta[i] = y;
 179                 linfo->Psi[i]   = z;
 180                 fprintf(stderr, "%f %f %f\n",linfo->Phi[i]*DEGREE, linfo->Theta[i]*DEGREE, linfo->Psi[i]*DEGREE);
 181         }
 182 }
 183
 184 void
 185 lProjectionDirectionMapCreate(lProjectionDirectionMapCreateInfo* linfo, int mode)
 186 {
 187         int    i, j, k;
 188         double radius;
 189         int flag;
 190
 191         linfo->x = MemoryAllocate(double, linfo->LineNum, "in lProjectionDirectionMapCreate");
 192         linfo->y = MemoryAllocate(double, linfo->LineNum, "in lProjectionDirectionMapCreate");
 193         j = 0;
 194         k = 0;
 195         for(i=0; i<linfo->LineNum; i++){
 196                 flag = 0;
 197                 switch(linfo->mode){
 198                         case 0 :
 199                                 radius  = linfo->Theta[i]*DEGREE;
 200                                 if(radius < 0)   radius += 360;
 201                                 if(radius > 180) radius -= 180;
 202                                 if(radius > 90)   radius = 180 - radius;
 203
 204                                 radius = 90 - radius;
 205
 206                                 linfo->x[i] = radius*cos((double)linfo->Phi[i]);
 207                                 linfo->y[i] = radius*sin((double)linfo->Phi[i]);
 208                         break;
 209                         case 1:
 210                                 if(linfo->Theta[i]*DEGREE < 0) radius = linfo->Theta[i]*DEGREE + 360;
 211                                 else                           radius = linfo->Theta[i]*DEGREE;
 212
 213                                 if(radius > 180) radius -= 180;
 214                                 if(radius > 90){
 215                                         radius = 180 - radius;
 216                                         flag = 1;
 217                                 }
 218
 219                                 radius = 90 - radius;
 220
 221                                 if(flag == 0){
 222                                         linfo->x[i] = radius*cos((double)(linfo->Phi[i]/2));
 223                                         linfo->y[i] = radius*sin((double)(linfo->Phi[i]/2));
 224                                 }
 225                                 else{
 226                                         linfo->x[i] = radius*cos((double)(linfo->Phi[i]/2+M_PI));
 227                                         linfo->y[i] = radius*sin((double)(linfo->Phi[i]/2+M_PI));
 228                                 }
 229                         break;
 230                         case 2 :
 231                                 linfo->x2 = MemoryAllocate(double, linfo->LineNum, "in lProjectionDirectionMapCreate");
 232                                 linfo->y2 = MemoryAllocate(double, linfo->LineNum, "in lProjectionDirectionMapCreate");
 233                                 if(linfo->Theta[i]*DEGREE < 0) radius = linfo->Theta[i]*DEGREE + 360;
 234                                 else                           radius = linfo->Theta[i]*DEGREE;
 235
 236                                 if(radius > 180) radius -= 180;
 237                                 if(radius > 90){
 238                                         radius = 180 - radius;
 239                                         flag = 1;
 240                                 }
 241
 242                                 radius = 90 - radius;
 243
 244                                 if(flag == 0){
 245                                         linfo->x[j] = radius*cos((double)(linfo->Phi[i]));
 246                                         linfo->y[j] = radius*sin((double)(linfo->Phi[i]));
 247                                         j++;
 248                                 }
 249                                 else{
 250                                         linfo->x2[k] = radius*cos((double)(linfo->Phi[i]));
 251                                         linfo->y2[k] = radius*sin((double)(linfo->Phi[i]));
 252                                         k++;
 253                                 }
 254                         break;
 255                         case 3 :
 256                         break;
 257                 }
 258         }
 259         switch(linfo->mode){
 260                 case 2 :
 261                         linfo->numP = j;
 262                         linfo->numM = k;
 263                 break;
 264         }
 265 }
 266
 267 void
 268 __OutputFileWrite(lProjectionDirectionMapCreateInfo* linfo, FILE* fpt, int mode)
 269 {
 270         int i;
 271         int max;
 272
 273         switch(linfo->mode){
 274                 case 0 :
 275                 case 1 :
 276                         for(i=0; i<linfo->LineNum; i++){
 277                                 fprintf(fpt, "%6.2f,%6.2f\n", linfo->x[i], linfo->y[i]);
 278                         }
 279                 break;
 280                 case 2 :
 281                         if(linfo->numP > linfo->numM){
 282                                 max = linfo->numP;
 283                                 for(i=0; i<max; i++){
 284                                         if(i > linfo->numM){
 285                                                 fprintf(fpt, "%6.2f,%6.2f\n",linfo->x[i], linfo->y[i]);
 286                                         }
 287                                         else{
 288                                                 fprintf(fpt, "%6.2f,%6.2f \t %6.2f,%6.2f\n",linfo->x[i], linfo->y[i], linfo->x2[i], linfo->y2[i]);
 289                                         }
 290                                 }
 291                         }
 292                         else{
 293                                 max = linfo->numM;
 294                                 for(i=0; i<max; i++){
 295                                         if(i > linfo->numP){
 296                                                 fprintf(fpt, "              \t %6.2f,%6.2f\n",linfo->x2[i], linfo->y2[i]);
 297                                         }
 298                                         else{
 299                                                 fprintf(fpt, "%6.2f,%6.2f \t %6.2f,%6.2f\n",linfo->x[i], linfo->y[i], linfo->x2[i], linfo->y2[i]);
 300                                         }
 301                                 }
 302                         }
 303                 break;
 304                 case 3 :
 305                         for(i=0; i<linfo->LineNum; i++){
 306                                 fprintf(fpt, "%7.2f,%7.2f\n", linfo->x[i], linfo->y[i]);
 307                         }
 308                 break;
 309         }
 310 }