4t1c.cc

// 4t1c is C++ generator of curves at the image http://mizugadro.mydns.jp/t/index.php/File:4t1.jpg
// 4t1c is C++ generator of curves at the image http://tori.ils.uec.ac.jp/TORI/index.php/File:4t1.jpg
// In order to combine these curves with data by Mainichi, the file 4t1.tex is also required
// Copyleft 2011 by Dmitrii Kouznetsov.
// At the reuse, please, indicate the source; this helps to trace errors, of any. 

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#define DB double
#define DO(x,y) for(x=0;x<y;x++)
using namespace std;
#include <complex>
typedef complex<double> z_type;
#define Re(x) x.real()
#define Im(x) x.imag()
#define I z_type(0.,1.)
#include "conto.cin"
double A[23][4]=
 {{ -26.8, 125.0,   3.8, 1},
  {  -7.5,  81.4,   6.2, 2},
  { -39.0,  80.0,   6.2, 3},
  { -39,    33.4,   9.8, 4},
  { -11,    7.8, 1770.7, 5},
  { -61.5,  7.5,    8.5, 6},
  {  -3.5, -16.2, 299.7, 7},
  { -28.5, -34.5,  34.1, 8},
  {-110.3, -39.5,   8.8, 9},
  {  -5,   -42.5,  33.2, 10},
  { -80,   -45,    11.8, 11},
  { -49.5, -46,    17.4, 12},
  {-125,   -51.5,   8.7, 13},
  { -34,   -61.2,  15.1, 14},
  { -31,   -68.8,  10.8, 15},
  {-132.5, -71,     8.3, 16},
  { -15,   -73,     8.4, 17},
  { -60.8, -74.3,   7.9, 18},
  { -31.9, -82.9,  10.2, 19},
  {-107.5, -86.1,  10.8, 20},
  {-107.6, -97.6,   7.1, 21},
  { -64.6,-100.6,   7.0, 22}
  }; // Length is measured in triks; roughly, 1trik=3km
DB T[23],U[23],F[23],W[23];
DB g(DB u, DB v, DB a, DB b){u/=a; v/=b; return  1./cosh(sqrt(u*u+v*v));}

DB X0=63.; DB Y0=90.;
DB g5(DB x, DB y, DB t) { x-=X0; y-=Y0; DB c=cos(t); DB s=sin(t);
                       return g(c*x+s*y, -s*x+c*y, 27.,3.9); }

DB g4(DB x, DB y, DB t) { DB c=cos(t); DB s=sin(t);
                       return g(c*x+s*y, -s*x+c*y, 17.,3.9); }

DB g45(DB x, DB y){ return 50000.*g5(x,y,-2.341)+1100.*g4(x,y,-1.98);}

main(){ int j,k,m,n; DB x,y,z, p,q, t;
//DO(k,23) { T[k]=A[k][0]; U[k]=A[k][1]; F[k]=A[k][2];
//     printf("%2d %6.1f %6.1f %6.1f\n",k,T[k],U[k],F[k]);}

int M=300,M1=M+1;
int N=400,N1=N+1;
DB X[M1],Y[N1], g[M1*N1], w[M1*N1]; // w is working array.
char v[M1*N1]; // v is also working array

FILE *o;o=fopen("4t1c.eps","w");ado(o,250,250);
fprintf(o,"140 110 translate\n");
DO(m,M1) X[m]=-200+1.*m;
DO(n,N1) Y[n]=-110+1.*n;
// for(m=-100;m<101;m+=100){    M(m, -100)L(m,100)}
// for(n=-100;n<201;n+=100){    M(-100,n)L(100.,n)}
// fprintf(o,".006 W 0 0 0 RGB S\n"); //frame centered at Fukushima 1 plant

DO(m,M1){x=X[m]; //printf("run at x=%6.3f\n",x);
DO(n,N1){y=Y[n]; g[m*N1+n]=g45(x,y); }}

fprintf(o,"1 setlinejoin 1 setlinecap\n");  q=800.;
//conto(o,g,w,v,X,Y,M,N, 4. ,-q,q); fprintf(o,".1 W 0 0 0 RGB S\n");
conto(o,g,w,v,X,Y,M,N, 10. ,-q,q); fprintf(o,"3 W 0 0 1 RGB S\n");
//conto(o,g,w,v,X,Y,M,N, 30. ,-q,q); fprintf(o,".1 W 0 0 0 RGB S\n");
conto(o,g,w,v,X,Y,M,N, 100. ,-q,q); fprintf(o,"3 W 0 .9 0 RGB S\n");
//conto(o,g,w,v,X,Y,M,N, 300. ,-q,q); fprintf(o,".1 W 0 0 0 RGB S\n");
conto(o,g,w,v,X,Y,M,N,1000. ,-q,q); fprintf(o,"3 W 1 0 0 RGB S\n");
fprintf(o,"showpage\n%cTrailer",'%'); fclose(o);
system("epstopdf 4t1c.eps");
system(    "open 4t1c.pdf");
}

// //