% script to generate from X-Y-W(X,Y)-Z(X,Y)-Q(X,Y) data regular grids of 
%                         X-Z-W(X,Z)-Y(X,Z)-Q(X,Z) type, currently only
% standard Legendre transform permutatations are generated, see eof for 
% list of tables.
%
%                                   JADC, 10/2008
%
% unit conversion on output assumed for Y, bar->Pa at end of file otherwise no unit 
% assumptions. 
%                        
%                                 interpolation method
iplot = 2;
jplot = 4;
ncont = 20;
%method = 'cubic'; % is unstable for STP table from PTVS input
method = 'linear';
%method = 'v4'; % good but too much memory
%method = 'nearest'; % no good

% read x-y-z array where z is the conjugate to y 
data_file=uigetfile('c*pl*', 'select y-conjugate x-y-z file (werami grid format)');
enter = fopen(data_file, 'rt');filename = fscanf(enter, '%s', 1);nx = fscanf(enter, '%f', 1);ny = fscanf(enter, '%f', 1);xmin = fscanf(enter, '%f', 1);ymin = fscanf(enter, '%f', 1);xinc = fscanf(enter, '%f', 1);yinc = fscanf(enter, '%f', 1);xname = fscanf(enter, '%s', 1);yname = fscanf(enter, '%s', 1);fclose(enter);
Z = rot90 (fliplr (reshape(textread(data_file, '%f','headerlines',4),nx,ny)));
% read x-y-assemblage array
data_file=uigetfile('c*pl*', 'select x-y-assemblage index file (werami grid format)');
enter = fopen(data_file, 'rt');filename = fscanf(enter, '%s', 1);nx = fscanf(enter, '%f', 1);ny = fscanf(enter, '%f', 1);xmin = fscanf(enter, '%f', 1);ymin = fscanf(enter, '%f', 1);xinc = fscanf(enter, '%f', 1);yinc = fscanf(enter, '%f', 1);xname = fscanf(enter, '%s', 1);yname = fscanf(enter, '%s', 1);fclose(enter);
P = rot90 (fliplr (reshape(textread(data_file, '%f','headerlines',4),nx,ny)));
% read x-y-w array where w is the conjugate to x
data_file=uigetfile('c*pl*', 'select x-conjugate x-y-w file (werami grid format)');
enter = fopen(data_file, 'rt');filename = fscanf(enter, '%s', 1);nx = fscanf(enter, '%f', 1);ny = fscanf(enter, '%f', 1);xmin = fscanf(enter, '%f', 1);ymin = fscanf(enter, '%f', 1);xinc = fscanf(enter, '%f', 1);yinc = fscanf(enter, '%f', 1);xname = fscanf(enter, '%s', 1);yname = fscanf(enter, '%s', 1);fclose(enter);
W = rot90 (fliplr (reshape(textread(data_file, '%f','headerlines',4),nx,ny)));
%                                 make regularly spaced X, Y, Z coordinate vectors
mnz = min(Z(:));mxz = max(Z(:));dz = (mxz - mnz)/(ny-1); z = mnz:dz:mxz;
mnw = min(W(:));mxw = max(W(:));dw = (mxw - mnw)/(ny-1); w = mnw:dw:mxw;
x = xmin:xinc:xmin+(nx-1)*xinc; y = ymin:yinc:ymin+(ny-1)*yinc; 
%                                  
%                                 make 2 grid for the X-Y coordinates of at 
%                                 which Z is known and a regularly spaced X-Z, W-Y
%                                 and W-Z grids
[X,Y] = meshgrid(x,y);
[XI,ZI] = meshgrid(x,z);
[WII,YII] = meshgrid(w,y);
[WIII,ZIII] = meshgrid(w,z);
%                                 __________________________________
%                                 T-V coordinate frame
%                                 the interpolated value of Y is:
YI = griddata(X(:),Z(:),Y(:),XI,ZI,method); 
subplot (iplot,jplot,1)
contourf(XI,ZI,YI,ncont);colorbar;title('T-V-pressure')
hold on
%plot3(X,Z,Y,'o', 'MarkerSize',1)
%                                 phase relations
PI = griddata(X(:),Z(:),P(:),XI,ZI,method); 
%subplot (iplot,jplot,2)
%contourf(XI,ZI,PI,ncont);colorbar;title('T-V-phase')
%                                 __________________________________
%                                 S-P coordinate frame
%                                 the interpolated value of X is:
XII = griddata(W(:),Y(:),X(:),WII,YII,method); 
ZII = griddata(W(:),Y(:),Z(:),WII,YII,method); 
subplot (iplot,jplot,2)
contourf(WII,YII,XII,ncont);colorbar;title('S-P-temperature')
hold on
%plot3(W,Y,X,'o', 'MarkerSize',1)
%                                 phase relations
PII = griddata(X(:),Z(:),P(:),XI,ZI,method); 
%subplot (iplot,jplot,4)
%contourf(WII,YII,PII,ncont);colorbar;title('S-P-phase')
%                                 __________________________________
%                                 S-V coordinate frame
%                                 the interpolated value of X & Y are:
XIII = griddata(W(:),Z(:),X(:),WIII,ZIII,method); 
YIII = griddata(W(:),Z(:),Y(:),WIII,ZIII,method);
PIII = griddata(W(:),Z(:),P(:),WIII,ZIII,method);
%subplot (iplot,jplot,5)
%contourf(WIII,ZIII,PIII,ncont);colorbar;title('S-V-phase')
subplot (iplot,jplot,3)
contourf(WIII,ZIII,YIII,ncont);colorbar;title('S-V-pressure')
hold on
%plot3(W,Z,Y,'o', 'MarkerSize',1)
subplot (iplot,jplot,4)
contourf(WIII,ZIII,XIII,ncont);colorbar;title('S-V-temperature')
hold on
%plot3(W,Z,X,'o', 'MarkerSize',1)
%                                 _________________________________
%                                 raw data
subplot (iplot,jplot,5)
contourf(X,Y,P,ncont);colorbar;title('P-T-phase')
subplot (iplot,jplot,6)
contourf(X,Y,W,ncont);colorbar;title('P-T-S')
subplot (iplot,jplot,7)
contourf(X,Y,Z,ncont);colorbar;title('P-T-V')

TTP = X; PTP = Y*1d5; STP = W; VTP = Z; QTP = P;
TTV = XI; VTV = ZI; PTV = YI*1d5; QTV = PI;
SSP = WII; PSP = YII*1d5; TSP = XII; VSP = ZII; QSP = PII;
SSV = WIII; VSV = ZIII; TSV = XIII; PSV = YIII*1d5; QSV = PIII;

save table TTP PTP STP VTP QTP TTV VTV PTV QTV SSP PSP TSP VSP QSP SSV VSV TSV PSV QSV