###########################################################################################
#
#
#  OCTAVE SCRIPT TO:
#       - READ ALEX DATA FROM TEKTRONIK SCOPES,
#	- READ DATA FROM STREAK IMAGES,
#	- OBTAIN FROM THEM RADIUS, AND VELOCITY DATA AND
#	- TRANSFORM THESE DATA INTO ELECTRICAL ENERGY DATA AND RADIUS AND VELOCITY TO CHECK THEM.
#    Made by Gonzalo Rodríguez Prieto
#       (gonzalo#rprietoATuclm#es)
#       (Mail: change "#" by "." and "AT" by "@")
#              Version 2.45
#
#
#########################################################
#
#  It uses the functions: 
#                        chan
#                        supsmu
#                        baseline
#			 ratio
#			 deri
#			 poladj
#			 vreduce
#  They must be in the same directory.
#
###########################################################################################

more off; %To make the lines display when they appear in the script, not at the end of it.

clear; %Just in case there is some data in memory.

tic; %Total time of the script.


#########################################################################################################
#
#	MAIN PARAMETERS FOR ALL THE SCRIPT:
#
#########################################################################################################

#Streak camera sweep range (For later calibration)
sweep = 097.66; %nanosec/px (For 50 microseconds in 512 rows in the picture.) [50 000 / 512]

#Metal and wire parameters:
rho = 19250; %Kg/m³ Tungsten density.
lengthwire = 2.*1e-2; %m Wire length
r0 = 0.05e-3; %m (Wire radius)
B = 2; %Constant between 2 and 4 for the EOS
Bi = 411e6; %Pascals (Tungten bulk modulus)
mass = rho * lengthwire * pi * r0 * r0 %Kilograms Wire mass

#########################################################################################################
#
#	SCOPE DATA INTERPRETATION PART.
#
#########################################################################################################



###
# Reading information over the shot features:
###

#Choosing the shot:
disp("On the ALEX shoot name, use alway capital letters.")
#shot = input("Which ALEX shoot to transform?\n","s");
shot = "ALEX087"; %TO make it unique

#Reading the TXT file with info over the shot:
#Shot name transformation:
infoshot = horzcat(shot,".txt");
#Open it:
[fileinfo, msg] = fopen(infoshot, "r");
if (fileinfo == -1) 
   error ("Script alex: Unable to open file name: %s, %s",infoshot, msg); 
endif;
#Putting the information into a cell structure:
txtinfo = textscan(fileinfo,"%s", "treatasempty","/");
fclose(fileinfo); #Closing the file
txtinfo = cell2mat(txtinfo); #To convert the cell structure in a matrix of strings.

#Putting the channels info in matrices of strings:
#(This is valid if the format of file SHOT.txt follows the procedure of ALEX data)
channels = ""; #Channel name
chantype = ""; #Channel probe type
for i=1 : length(txtinfo)
 if (cell2mat(strfind(txtinfo(i),"CH")) ==1) %We have the channels here, with the type following.
   % Standard format of the *.txt file for ALEX.
   channels = strvcat(channels,txtinfo(i));
   chantype = strvcat(chantype,txtinfo(i+1));
 elseif(cell2mat(strfind(txtinfo(i),"Del")) ==1) %We should stop to avoid channels confusion.
   #After channels info cames other shot info, not stored or used in this version of ALEX.m
   break;
 endif;
endfor;


###
#Taking all channel data in one matrix (VOL): 
###

for i=1 : rows(channels)
 channame = strcat(shot,".",channels(i,:),".CSV");
 fprintf(stdout,"Taking data from %s\n", channame);
 [t(:,i),vol(:,i),points(i)] = chan(channame); #chan is the function to take the channel file CHANNAME info.
endfor;


###
#Transforming scope data from (t,volts) into (t,dimensional units) taking into account the type:
###

#Selecting the variable type from the file information:
selection = []; #Defining the variable
for i=1 : rows(channels)
#Choosing which adaptation is necessary from the channel type characters:
if strcmp(chantype(i,1:4),"Div.")==1
 selection = [selection,1]; %Number 1: Resistive divider
endif;
if strcmp(chantype(i,1:5),"Sonda")==1
 selection = [selection,2]; %Number 2: Tektronik voltage probe
endif;
if strcmp(chantype(i,1:3),"Rog")==1
 selection = [selection,3]; %Number 3: Rogowsky coil probe
endif;
if ( length(chantype(i,:))>8 )
    if strcmp(chantype(i,1:10),"Photodiode")==1
     selection = [selection,4]; %Number 4: Photodiode signal
    endif;
endif;
endfor;

# Data transformation of the channels from previous selection:
for i=1 : columns(selection)
 switch(selection(i))
    case 1 %My resistive divider data
      volt(:,2) = vol(:,i) ./ 0.00032; %Wire volts
      volt(:,1) = t(:,i).*1e6; %time in us
      volt(:,2) = supsmu(volt(:,1),volt(:,2),"span",0.01); %Smoothing the voltage data.
	#It allows for a better treatment of the data, without losing too much information.
      disp("Transforming voltage data from resistive divider");
    case 2 %The tektronik probe
      volt(:,2) = vol(:,i); %Wire volts. No transformation, but not clear if correct.
      volt(:,1) = t(:,i).*1e6; %time in us
      volt(:,2) = supsmu(volt(:,1),volt(:,2),"span",0.01); %Smoothing
      disp("Transforming voltage data from Tektronik probe");
    case 3 %Rogowsky coil
      dint(:,2) = (vol(:,i).*63.095)./5.85e-9; 
        %Following calibration page 33 of "Logbook of diverse projects 01"
      dint(:,1) = t(:,i).*1e6; %time in us
      dint(:,2) = supsmu(dint(:,1),dint(:,2),"span",0.01); %Smoothing
      disp("Transforming Rogowsky signal");
    case 4 %Photodiode signal (Normalized to 1)
      phot(:,2) = vol(:,i);
      phot(:,2) = phot(:,2)./max(phot(:,2)); 
      phot(:,1) = t(:,i).*1e6; %time in us
      phot(:,2) = supsmu(phot(:,1),phot(:,2),"span",0.02); %Smoothing
      disp("Transforming photodiode signal");
 endswitch
endfor;

#Checking and warning for missed data on the shot:
if (exist("volt","var")==0)
  warning("Script alex: There is no voltage signal in shot %s",shot);
elseif (exist("dint","var")==0)
  warning("Script alex: There is no intensity signal in shot %s",shot);
elseif (exist("phot","var")==0)
  warning("Script alex: There is no photodiode signal in shot %s",shot);
endif;


###
#Removing values on the baseline for integration or further operations of the data:
###

disp("Removing the baselines form the scope data.");
#Finding the baseline values (if the variables exist):
#Putting the voltage baseline to zero:
if (exist("volt","var")==1)
  [basev,mediav,puntosv] = baseline(volt(1:250,2)); %Voltage signal
  if (puntosv<=10)
    warning("Script alex: only %u points in voltage baseline",puntosv);
  endif;
  volt(:,2) = volt(:,2) - mediav;
endif;
#The Rogowsky signal baseline to zero:
if (exist("dint","var")==1)
  [basei,mediai,puntosi] = baseline(dint(1:250,2)); %Intensity signal
  if (puntosi<=10)
    warning("Script alex: only %u points in intensity derivative baseline",puntosi);
  endif;
  dint(:,2) = dint(:,2) - mediai;
endif;
#Photodiode signal to zero:
if (exist("phot","var")==1)
  [basep,mediap,puntosp] = baseline(phot(1:200,2)); %Photodiode signal
  if (puntosp<=10)
    warning("Script te: only %u points in photodiode baseline",puntosp);
  endif;
  phot(:,2) = phot(:,2) - mediap;
endif;


###
# Integrating the intensity derivative:
###

if (exist("dint","var")==1)
  int = - cumtrapz(dint(:,1)./1e6,dint(:,2)); %The minus because of the signal in the voltage and the derivative.
  #The division of 1e6 is because we need to use ALL unit in I.S!!!!!
endif;


####
## Finding the electrical power and energy delivered by time and the total: 
####

if ( (exist("int","var")==1) && (exist("volt","var")==1) ) 
  elpowr = volt(:,2).*int; #Electrical power (Watts)
  #Electrical energy delivered:
  elenergia(:,1) = dint(:,1); %Time (microseconds)
  elenergia(:,2) = cumtrapz(elpowr).*(abs(dint(1,1)-dint(2,1)).*1e-6); %Energy (Joules)
endif;







#########################################################################################################
#
#	STREAK IMAGE INTERPRETATION PART.
#
#########################################################################################################


####
## Opening the picture and choosing some parameters.
####

filename = horzcat(shot,".dat"); %Compatibility with previous code and data treatment.

disp("Opening streak image.");

[file, msg] = fopen(filename, "r");
if (file == -1) 
   error ("Unable to open file name: %s, %s",filename, msg); 
endif; 

#Standard deviation parameter (Value to make the binarization)
std_dev_par = 20; %This seems to work, change if results are very bad.

#Charge the image as a matrix: (IT MUST BE ALREADY FORMATTED AS SUCH)
T = dlmread(file);
fclose(file);

#The binarized matrix. empty with zeroes now:
Out = zeros(rows(T),columns(T));

#The output vector. As before, it starts with zeroes inside:
move = zeros(rows(T),8);

#Maximum of the matrix/2. For showing the results in images and finding edges.
valout = max(max(T))/2;


###
# Making the binarized matrix.
###
#   Loop to fill the OUT matrix with the positions and binarized image: 
#  It places "valout" on the output matrix. So this matrix has only two values: 0
# and "valout". Then it is easy to look for the external edges.
#  It is done on a double loop because with vector logics and so on,
# the program gives not right results.

disp("Making the binarized matrix.");

for i = 1 : rows(T) #For every row of the picture...
    row = T(i,:); # put the row on a variable.
    if ( std(row) > std_dev_par ) #When there is data on the row, no just noise ->
    % Signaled by a standard deviation higher than "std_dev_par".
      halfrow = max(row)/2; # Take half of the maximum value for this row and use it as mark for margins
	for j = 1 : columns(T)
           %  In every column of the OUT row change the value to "valout" if 
           % the column value is higher than halfrow.
	   if ( row(j) > halfrow )
	     Out(i,j) = valout;
	   endif;
	endfor;
    endif;
endfor;


###
# Find the positions of the edges in the binary matrix
###
#Loop for finding the real positions pixels table:

disp("Finding the edges in the matrix.");

j = 0; %Initialize the "j" variable.
for i = 1 : rows(Out) %For every row of the picture...
   pos = find(Out(i,:)==valout); %Find the positions of "valout" in the output matrix and
   # store them in an index vector(POS)
   if (columns(pos) > 1) %When there are really values
     j = j+1;
     move(j,1) = i; %First column is time (In pixels units)
     move(j,2) = pos(1); %Second column is space on "first" side
     move(j,3) = pos(end); %Third column is space on the "last" side
     if (j==1) %Finding the center position for the first column:
       zenter = abs(move(1,2)-move(1,3))/0.5 + min(move(1,2),move(1,3));
     endif;
   endif;
endfor;


###
# Streak data conversion in radius and velocity over time.
###

disp("Streak data treatment and transformation into radius and velocity over time.");

#Remove the zeros on the matrix results:
 move = move( (move(:,2)!=0),:); %Take out the positions were the second column has zero value.
#(Made using logical indexing)

#Calibration of pixels in time and space and centering:
#CAREFUL: THIS CALIBRATION IS VALID FROM ALEX086 TO ALEX099 ONLY !!!!!
move(:,4) = move(:,1) .* sweep; #Passing pixels to nanoseconds.
move(:,5) = (move(:,2)-zenter) .* 0313; #Passing pixels to micrometers and centering.
move(:,6) = (move(:,3)-zenter) .* 0313; #Passing pixels to micrometers and centering.

#smoothing radius data with the function "supsmu", check function help to see how it works:
move(:,5) = supsmu(move(:,4),move(:,5),"span",0.01);
move(:,6) = supsmu(move(:,4),move(:,6),"span",0.01);


###
# Data transforming into radius expansion and its velocity
###

#Putting in zero the radius displacement (different for each side)
#This configuration depends on the streak image orientation. 
#Is now made for "up" left side and "down" right side.
move(:,5) = -(move(:,5) - max(move(:,5))); 
move(:,6) = move(:,6) - min(move(:,6)); 

#Name the displacement vector in a more convinient way:
rup = ( (move(:,5) .* 1e-6) + r0) ./ r0; %transforming micrometers into dimensionless space.
rdown = ( (move(:,6) .* 1e-6) + r0) ./ r0; %transforming micrometers into dimensionless space.
t = move(:,4) .* 1e-9; %transforming nanoseconds in seconds.

#Deriving to obtain velocity (in micrometers/nanosecond):
  h = abs(move(10,4)-move(11,4));
dev = deri(move(:,5),h);
dev2 = deri(move(:,6),h);

move(1:columns(dev),7) = 1000 .* dev; %Transforming it in m/s from um/ns and placing it in the "move" matrix.
move(1:columns(dev),8) = 1000 .* dev2; %Transforming it in m/s from um/ns and placing it in the "move" matrix.
#Renaming velocities:
vup = move(:,7); 
vdown = move(:,8);

#Remove zeroes:
tv = t(vup!=0);
vup(vup==0) = [];
vdown(vdown==0) = [];

#Adjusting to a polynomium the radius to derive it to find the velocity:
[ruppol,error_ruppol,rupol_estr] = polyadj(t,rup.*r0,0.90); %L[t] meters(s) Radius as a polynomial.
vup2 = polyval(polyder(ruppol),t); %m/s. Evaluating velocity as derivative of the polinomial radius over time.



plot(tv,vup,"*r","markersize",10,t,vup2,"-k")
title(shot);
legend({"Velocity as derivative","velocity from polynomial radius"});
ylabel("m/s");
xlabel("t (mseconds)");
print(horzcat(shot,"-velocity.jpg"),"-djpg");


plot(t.*1e6,rup.*r0.*1e3,"*k",t.*1e6,polyval(ruppol,t).*1e3,"--r")
axis([0 1 0  (max(rup.*r0.*1e3)+0.5)]); %Axis limits. First x axis, then y axis (In the units of the vector files).
axis("autox");
title(shot);
legend({"Radial points","Polinomialfitting"});
ylabel("Radius (mm)");
xlabel("time (microseconds)");
print(horzcat(shot,"-radius.jpg"),"-djpg");

%Writing down the binarized streak: 
imwrite(Out,"ALEX087-bin-streak.jpg","jpg")
%Writting down the streak:
figure;
imagesc(T); %Scaling image form the data matrix.
print(horzcat(shot,"-Streak.jpg"),"-djpg");


#####
## Plotting second figure (Voltage, Intensity).
#####
disp("Plotting figure 2...");
intensity = [dint(:,1),int];

#Reducing number of points. Too many points make the Octave slow.
	volt = vreduce(200,volt);
	intensity = vreduce(200,intensity);

lengthax = ( max(intensity(:,2)) + 500) .*1e-3; %Dynamical axis length.
 
[ax hint hvolt] =plotyy(intensity(:,1),intensity(:,2).*1e-3, %Intensity data
volt(:,1),volt(:,2).*1e-3); %Voltage data.
axes(ax(1)); %Intensity axis
set(hint,"linestyle",".-","markersize",3,"colorline","red");
ylabel("Intensity(kA)");
axes(ax(2)); %Voltage axis
set(hvolt,"linestyle","--","markersize",7,"colorline","blue");
ylabe("Voltage(kV");

text(6.5,round(lengthax/2),'.- intensity (kA)','color','red');  %Labels on the plot. First options, coordinates on the plot dimensions
text(6.5,round(lengthax/2.25),'-- voltage (kV)','color','blue');
xlabel("time (microseconds)");
title(horzcat(shot," Electric data."));
axis([-2 10 -lengthax lengthax]); %Axis limits. First x axis, then y axis (In the units of the vector files).
print(horzcat(shot,"-electric.jpg"),"-djpg");

[rdownpol,error_rdownpol,rdownpol_estr] = polyadj(t,rdown.*r0,0.90); %L[t] meters(s) Radius as a polynomial.
vdown = polyval(polyder(rdownpol),t); %m/s Velocity as derivative of radius.


more on; #Revert more normal comments behaviour.

###
# Total processing time
###
timing = toc;
disp("Script alex execution time:")
disp(timing)
disp(" seconds")

#That's...that's all, folks! 
