clf;
clear;

file = 'shape1';
fid  = fopen(file);

mesh_t=50;
mesh_phi=50;
maxpoly=10;

axmin = fscanf(fid,'%f',1);
axmax = fscanf(fid,'%f',1);
aymin = fscanf(fid,'%f',1);
aymax = fscanf(fid,'%f',1);
azmin = fscanf(fid,'%f',1);
azmax = fscanf(fid,'%f',1);
maxs = fscanf(fid,'%d',1);
axis([axmin axmax aymin aymax azmin azmax])
axis('equal');
axis('image');

xlabel('X');
ylabel('-Z');
zlabel('Y');
hold on;

axis off;

view([-0.5 -1.2 0.7]);

t = (0:mesh_t)/mesh_t;
tn(1,:) = 0.d0 * ones(1,mesh_t+1);
tn(2,:) = 1.d0 * ones(1,mesh_t+1);
for n=3:maxpoly+1
 tn(n,:) = tn(n-1,:).*(t-0.5);
end

phi = (0:mesh_phi)/mesh_phi*2*pi;
sinn = sin(phi);
cosn = cos(phi);

for is=1:maxs
  mpoly_x = fscanf(fid,'%d',1);
  ax = fscanf(fid,'%f',mpoly_x);
  mpoly_y = fscanf(fid,'%d',1);
  ay = fscanf(fid,'%f',mpoly_y);
  mpoly_z = fscanf(fid,'%d',1);
  az = fscanf(fid,'%f',mpoly_z);
  mpoly_r = fscanf(fid,'%d',1);
  ar = fscanf(fid,'%f',mpoly_r);

  coeff = 1:mpoly_x;
  xx = ax' * tn(coeff+1,:);
  cx = (coeff-1) .* ax' * tn(coeff,:);
  coeff = 1:mpoly_y;
  yy = ay' * tn(coeff+1,:);
  cy = (coeff-1) .* ay' * tn(coeff,:);
  coeff = 1:mpoly_z;
  zz = az' * tn(coeff+1,:);
  cz = (coeff-1) .* az' * tn(coeff,:);
  coeff = 1:mpoly_r;
  radius = ar' * tn(coeff+1,:);

  c2 = cx.*cx+cy.*cy+cz.*cz;
  c1 = sqrt(c2);
  cx = cx./c1;
  cy = cy./c1;
  cz = cz./c1;

% generate random reference systems for all solenoids
  vx = rand;
  vy = rand;
  vz = rand;

% take a cross-product v*c
  aax = vy*cz-cy*vz;
  aay = cx*vz-vx*cz;
  aaz = vx*cy-cx*vy;
  aa1 = sqrt(aax.*aax+aay.*aay+aaz.*aaz);
  aax = aax ./ aa1 .* radius;
  aay = aay ./ aa1 .* radius;
  aaz = aaz ./ aa1 .* radius;
% take the cross-product c*aa
  bbx = cy.*aaz-aay.*cz;
  bby = aax.*cz-cx.*aaz;
  bbz = cx.*aay-aax.*cy;

  xs = ones(mesh_phi+1,1)*xx + ( cosn' * aax + sinn' * bbx ); 
  ys = ones(mesh_phi+1,1)*yy + ( cosn' * aay + sinn' * bby );
  zs = ones(mesh_phi+1,1)*zz + ( cosn' * aaz + sinn' * bbz );
 
  mesh(xs,-zs,ys);
end;

box(axmin,-azmin,aymin,axmax,-azmax,aymax,'-.')

fclose(fid); 


