clear;
format('e',12);

// *** PREMのテーブルを読み出し ***
X = fscanfMat('PREM.txt');
Radius = 1E3 * X(:,1);		// 半径 (m)
Vp = X(:,2);			// P波速度 (m/s)
Vs = X(:,3);			// S波速度 (m/s)
RHO = 1E3 * X(:,4);		// 密度 (g/m^3)
Ks = 1E12 * X(:,5);		// 断熱体積弾性率(Pa)
Mu = 1E12 * X(:,6);		// 剛性率(Pa)
Nu = 1E12 * X(:,7);		// 
Pressure = 1E12 * X(:,8);	// 圧力(Pa)
Gravity = X(:,9);		// 重力加速度 (m/s^2)

// 半径と密度の外核の部分だけ取り出し
OCR = Radius(9:21);		// 外核の半径 Outer Core Radius
OCD = RHO(9:21);		// 外核の密度 Outer Core Density

// *** 内挿 ***
OCRp = linspace(OCR(1),OCR($));
OCDp = interp1(OCR,OCD,OCRp,'spline');

// *** プロット ***
// 外核全領域のプロット
scf(0);
plot(1E-3 * OCRp, 1E-3 * OCDp,'-r');
plot(1E-3 * Radius(9:21),1E-3 * RHO(9:21),'ob');
ylabel("Density (kg/s^3)");
xlabel("Radius (km)");

// 2300-2700kmの領域のプロット
scf(1);
//xsetech([0,0,0.95,0.95]);
// スプライン曲線による補間 
plot(1E-3 * OCRp, 1E-3 * interp1(OCR,OCD,OCRp,'spline'),'-r');
// 線形補間
plot(1E-3 * OCRp, 1E-3 * interp1(OCR,OCD,OCRp,'linear'),'--b');
// 最近接データの値による補間
plot(1E-3 * OCRp, 1E-3 * interp1(OCR,OCD,OCRp,'nearest'),'-.g');
plot(1E-3 * Radius(9:21),1E-3 * RHO(9:21),'ob');
zoom_rect([2300,11050,2700,11350]);
legend(['spline','linear','nearest']);
ylabel("Density (kg/s^3)");
xlabel("Radius (km)");