synthe.c

射线追踪程序

/*
 * This file is part of tomo3d
 *
 * Copyright (C) 2002, 2003, Sebastien Judenherc <[e-mail]>
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307
 * USA
 *
 */
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include "synthe.h"

extern float ANISO1(float x, float y, float z, float A, float P);
extern float ANISO2(float x, float y, float z, float A, float P);
extern float ANISO3(float x, float y, float z, float A, float P);
extern float ANISO4(float x, float y, float z, float A, float P);
extern float ANISO_G1(float x, float y, float z, float A, float P);
extern float ANISO_G2(float x, float y, float z, float A, float P);
extern float ANISO_G3(float x, float y, float z, float A, float P);

void initaniso();
void initaniso2(float az1, float dip1, float az2, float dip2);

float CB3D(float x, float y, float z, float A, float P)
{ float w=2.0*M_PI/P;
  return(1.0+A*(cos(w*x)*cos(w*y)*cos(w*z)));
}

float CBX(float x, float y, float z, float A, float P)
{ float w=2.0*M_PI/P;
  return(1.0+A*cos(w*y)*cos(w*z));
}

float CBY(float x, float y, float z, float A, float P)
{ float w=2.0*M_PI/P;
  return(1.0+A*cos(w*x+w*z));
}

float CBZ(float x, float y, float z, float A, float P)
{ float w=2.0*M_PI/P;
  return(1.0+A*cos(w*x+w*y));
}

float STEPX(float x, float y, float z, float A, float P)
{ if (y<P) return(1.-A);
  return(1.+A);
}

float STEPY(float x, float y, float z, float A, float P)
{ if (x<P) return(1.-A);
  return(1.+A);
}

float TOTO(float x, float y, float z, float A, float P)
{ if (z<100.0) return(1.0);
  if (z>180.0) return(1.0);
  if ((x>-40.)&&(x<40.)&&(y<40.)&&(y>-40.)) return(1.+A);
  return(1.0);
}

float HORS(float x, float y, float z, float A, float P)
{ if (z<200.0) return(1.0);
  if (z>250.0) return(1.0);
  if ((x>-25.)&&(x<25.)&&(y<25.)&&(y>-25.)) return(1.+A);
  return(1.0);
}

float CMPLX(float x, float y, float z, float A, float P)
{ float X1,Z1,X2,Z2,T,Z0;
  if (z<30.0) {
     if ((x>30.0)&&(x<40.0)) return(1.0+A);
     if ((y>0.0)&&(y<10.0)) return(1.0-A);
     if ((y>-100.0)&&(y<-90.0)) return(1.0-A);
  }
  if (z>200.0) return(1.0);
  if (z<60.0) return(1.0);
  Z1=60.0;  X1=50.0;
  Z2=200.0; X2=-100.0;
  T=(Z2-Z1)/(X2-X1);
  Z0=Z1-T*X1;
  if (z>Z0+y*T) return(1.0);

  Z1=60.0;  X1=20.0;
  Z2=200.0; X2=-130.0;
  T=(Z2-Z1)/(X2-X1);
  Z0=Z1-T*X1;
  if (z<Z0+y*T) return(1.0);
  return(1.0-A);
}

float TEST(float x, float y, float z, float A, float P)
{ 
  if (z>60.) return(1.0);
  if (z<30.) {
          if (sqrt((x*x)+(y*y))<3) {
                  return(0.92);
          } else {
                  return(1.);
          }
  }
  if (z<60.) {
          if (sqrt((x*x)+(y*y))<15) {
                  return(0.92);
          }
  }
  return(1.);
}

struct synthe_t *parseVfunc(char *s1, char *s2, char *s3)
{ struct synthe_t *tmp;
  tmp=(struct synthe_t *)malloc(sizeof(struct synthe_t));
  if (!sscanf(s2,"%f",&(tmp->P1))) return(NULL);
  tmp->P1/=100.;
  if (!sscanf(s3,"%f",&(tmp->P2))) return(NULL);
  if (!strcmp(s1,"TEST")) {
     tmp->fun=TEST;
     return(tmp);
  }
  if (!strcmp(s1,"CB3D")) {
     tmp->fun=CB3D;
     return(tmp);
  }
  if (!strcmp(s1,"CBX")) {
     tmp->fun=CBX;
     return(tmp);
  }
  if (!strcmp(s1,"CBY")) {
     tmp->fun=CBY;
     return(tmp);
  }
  if (!strcmp(s1,"CBZ")) {
     tmp->fun=CBZ;
     return(tmp);
  }
  if (!strcmp(s1,"STEPX")) {
     tmp->fun=STEPX;
     return(tmp);
  }
  if (!strcmp(s1,"STEPY")) {
     tmp->fun=STEPY;
     return(tmp);
  }
  if (!strcmp(s1,"TOTO")) {
     tmp->fun=TOTO;
     return(tmp);
  }
  if (!strcmp(s1,"HORS")) {
     tmp->fun=HORS;
     return(tmp);
  }
  if (!strcmp(s1,"CMPLX")) {
     tmp->fun=CMPLX;
     return(tmp);
  }
  if (!strcmp(s1,"ANISO1")) {
     tmp->fun=ANISO1;
     initaniso();
     return(tmp);
  }
  if (!strcmp(s1,"ANISO2")) {
     tmp->fun=ANISO2;
     initaniso();
     return(tmp);
  }
  if (!strcmp(s1,"ANISO3")) {
     tmp->fun=ANISO3;
     initaniso();
     return(tmp);
  }
  if (!strcmp(s1,"ANISO4")) {
     tmp->fun=ANISO4;
     initaniso();
     return(tmp);
  }
  if (!strncmp(s1,"ANISO_G",strlen("ANISO_G"))) {
     float val1;
     sscanf(s1+strlen("ANISO_G1_"),"%f",&val1);
     if (s1[7]=='1') tmp->fun=ANISO_G1;
     	else if (s1[7]=='2') tmp->fun=ANISO_G2;
     	     else if (s1[7]=='3') tmp->fun=ANISO_G3;
	          else return(NULL);
     initaniso2(56.0,val1,40.0,0.0);
     return(tmp);
  }
  return(NULL);
}

float compVsyn(struct synthe_t * syn, float x, float y, float z)
{ return(syn->fun(x,y,z,syn->P1,syn->P2)); }

void printfsynlist()
{ fprintf(stderr,"CB3D     <Amp> <Per>\n");
  fprintf(stderr,"CBX      <Amp> <Per>\n");
  fprintf(stderr,"CBY      <Amp> <Per>\n");
  fprintf(stderr,"CBZ      <Amp> <Per>\n");
  fprintf(stderr,"STEPX    <Amp> <Offset>\n");
  fprintf(stderr,"STEPY    <Amp> <Offset>\n");
  fprintf(stderr,"TEST     <Amp> <Offset>\n");
  fprintf(stderr,"HORS     <Amp> <Offset>\n");
  fprintf(stderr,"CMPLX    <Amp> <Offset>\n");
  fprintf(stderr,"ANISO    <Amp> <maxdep>\n");
  fprintf(stderr,"ANISO_G1_0   <Amp> <maxdep>\n");
}