안녕하세요 C언어 쌩초보입니다... 선배가 C++ 6.0 에서 짜놓은 코드로 visual studio 2015에서 해석을 돌리려 하는데
자꾸만 오류가 발생해서 일주일 째 골머리르 썩고 있습니다. 부디 지식 나눔을 해주시어 답을 내려주세요 ㅠㅠ
// AB_MAP.cpp : Defines the entry point for the console application.
//
#include "stdafx.h"
#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <string.h>
#include <fstream>
#include <math.h>
#include <time.h>
int num_lon, num_lat, num_M;
int control_b = 10;
double W,E,S,N,dD;
double Ml,Mu,dM;
double aa[91][101];
double bb[91][101];
double probability(double[5], double, double, double, double, double*);
int return_p[4];
int iteration(int);
void main()
{
int i, j, m, ii;
int iter;
double del_b;
double lat, lon, b, RMSE;
double a10, v_old, v_new;
double RMSE_a, RMSE_b;
double* original1;
double* original2;
double* original3;
double* result1;
double* result2;
double* result3;
double* sm_b;
char temp[80];
return_p[1] = 500; return_p[2] =1000; return_p[3] = 2400;
W = 123; E = 131; S = 33; N=43;
dD = 0.1;
num_lon = int((E-W)/dD); num_lat = int((N-S)/dD);
Ml=4.5; Mu=7; dM=0.5;
num_M = int((Mu-Ml)/dM);
FILE* fin;
FILE* fin_o1;
FILE* fin_o2;
FILE* fin_o3;
FILE* fin_r1;
FILE* fin_r2;
FILE* fin_r3;
FILE* finput;
FILE* frmse;
FILE* fsmoothb;
///////////////////////////////////////////////////////////////////////
///////////////////////////////////////// Read Input data /////////////
fin = fopen("input.txt","r");
fscanf(fin,"%24c\n",temp);
fscanf(fin,"%d %lf %lf %lf\n",&iter,&del_b,&RMSE_a,&RMSE_b);
printf("Iteration Number = %d\n",iter);
///////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////
///////////////////////////////////////// Obtain Original Data ////////
sprintf(temp,"seismic_map_%d_original.dat",return_p[1]);
fin_o1 = fopen(temp,"r");
sprintf(temp,"seismic_map_%d_original.dat",return_p[2]);
fin_o2 = fopen(temp,"r");
sprintf(temp,"seismic_map_%d_original.dat",return_p[3]);
fin_o3 = fopen(temp,"r");
b = double(control_b);
if(control_b == 10)
{
fsmoothb = fopen("smooth_b.txt","r");
sm_b = new double[9001];
}
original1 = new double[9001];
original2 = new double[9001];
original3 = new double[9001];
for(j=1;j<=num_lat;j++)
{
for(i=1;i<=num_lon;i++)
{
m = (j-1)*num_lon + i;
fscanf(fin_o1,"%lf",&original1[m]);
fscanf(fin_o2,"%lf",&original2[m]);
fscanf(fin_o3,"%lf",&original3[m]);
if(control_b == 10) fscanf(fsmoothb,"%lf\n",&sm_b[m]);
}
}
fclose(fin_o1);
fclose(fin_o2);
fclose(fin_o3);
if(control_b == 10) fclose(fsmoothb);
///////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////
//////////////////////////////// Scale by Original / New data /////////
sprintf(temp,"RMSE_%.2f.txt",b);
frmse = fopen(temp,"w");
sprintf(temp,"copy ab_data_ceus_%.2f.dat ab_data_ceus.dat",b);
system(temp);
if(iter == 0)
{
sprintf(temp,"copy ab_data_ceus.dat ab_data_ceus_0.dat");
system(temp);
iteration(0);
sprintf(temp,"copy seismic_map_500.dat seismic_map_500_0.dat");
system(temp);
sprintf(temp,"copy seismic_map_1000.dat seismic_map_1000_0.dat");
system(temp);
sprintf(temp,"copy seismic_map_2400.dat seismic_map_2400_0.dat");
system(temp);
fin_r1 = fopen("seismic_map_500.dat","r");
fin_r2 = fopen("seismic_map_1000.dat","r");
fin_r3 = fopen("seismic_map_2400.dat","r");
result1 = new double[9001];
result2 = new double[9001];
result3 = new double[9001];
RMSE = 0;
for(j=1;j<=num_lat;j++)
{
for(i=1;i<=num_lon;i++)
{
m = (j-1)*num_lon + i;
fscanf(fin_r1,"%lf",&result1[m]);
fscanf(fin_r2,"%lf",&result2[m]);
fscanf(fin_r3,"%lf",&result3[m]);
RMSE = RMSE + fabs(original1[m]-result1[m])+fabs(original2[m]-result2[m])+fabs(original3[m]-result3[m]); // RMSE check
}
}
printf("RMSE_a = %lf\n",RMSE);
fprintf(frmse,"%.6f\n",RMSE);
fclose(fin_r1);
fclose(fin_r2);
fclose(fin_r3);
}
else
{
for(ii=1;ii<=iter;ii++)
{
if( ii != 1) system("copy ab_data_ceus_temp.dat ab_data_ceus.dat");
sprintf(temp,"copy ab_data_ceus.dat ab_data_ceus_%d.dat",ii);
system(temp);
iteration(ii); // Gain New Data
sprintf(temp,"copy seismic_map_500.dat seismic_map_500_%d.dat",ii);
system(temp);
sprintf(temp,"copy seismic_map_1000.dat seismic_map_1000_%d.dat",ii);
system(temp);
sprintf(temp,"copy seismic_map_2400.dat seismic_map_2400_%d.dat",ii);
system(temp);
RMSE = 0;
fin_r1 = fopen("seismic_map_500.dat","r");
fin_r2 = fopen("seismic_map_1000.dat","r");
fin_r3 = fopen("seismic_map_2400.dat","r");
result1 = new double[9001];
result2 = new double[9001];
result3 = new double[9001];
for(j=1;j<=num_lat;j++)
{
for(i=1;i<=num_lon;i++)
{
m = (j-1)*num_lon + i;
fscanf(fin_r1,"%lf",&result1[m]);
fscanf(fin_r2,"%lf",&result2[m]);
fscanf(fin_r3,"%lf",&result3[m]);
RMSE = RMSE + fabs(original1[m]-result1[m])+fabs(original2[m]-result2[m])+fabs(original3[m]-result3[m]); // RMSE check
}
}
printf("RMSE_a = %lf\n",RMSE);
fprintf(frmse,"%.6f\n",RMSE);
fclose(fin_r1);
fclose(fin_r2);
fclose(fin_r3);
if(RMSE < RMSE_a || ii == iter) break;
finput = fopen("ab_data_ceus_temp.dat","w");
for(j=1;j<=num_lat;j++)
{
lat = 33.0+(j-1)*dD;
for(i=1;i<=num_lon;i++)
{
lon = 123.0+(i-1)*dD;
m = (j-1)*num_lon + i;
if( aa[i][j] != 0 ) v_old = pow(10,log10(aa[i][j])-bb[i][j]*Ml);
else v_old = 0;
if(original1[m] < 0.001) a10 = 0;
else
{
if( v_old < 1e-10 ) a10 = 0;
else if( result1[m] == 0) a10 = aa[i][j];
else
{
if( original1[m] >= result1[m] && original2[m] >= result2[m] && original3[m] >= result3[m] )
v_new = v_old * original2[m] / result2[m];
else if( original1[m] < result1[m] && original2[m] < result2[m] && original3[m] < result3[m] )
v_new = v_old * original2[m] / result2[m];
else v_new = v_old;
if( original1[m] >= result1[m] && original3[m] < result3[m] )
sm_b[m] = sm_b[m] + del_b;
else if( original1[m] < result1[m] && original3[m] >= result3[m] )
sm_b[m] = sm_b[m] - del_b;
if( sm_b[m] >= 2.0) sm_b[m] = 2.0;
else if( sm_b[m] <= 0.2 ) sm_b[m] = 0.2;
if(control_b == 10) a10 = pow(10,log10(v_new)+sm_b[m]*Ml);
else a10 = pow(10,log10(v_new)+b*Ml);
}
}
if(a10 == 0)
{
if(control_b == 10)
fprintf(finput,"%.2f %.2f %.8f %.8f 0.00\n",lon,lat,a10,sm_b[m]);
else
fprintf(finput,"%.2f %.2f %.8f %.8f 0.00\n",lon,lat,a10,b);
}
else
{
if(control_b == 10)
fprintf(finput,"%.2f %.2f %.8f %.8f %.2f\n",lon,lat,a10,sm_b[m],Mu);
else
fprintf(finput,"%.2f %.2f %.8f %.8f %.2f\n",lon,lat,a10,b,Mu);
}
}
}
fclose(finput);
}
}
fclose(frmse);
//////////////////////////////////////////////////////////////////
}
int iteration(int iter)
{
int i,ii,j,k1,k2,k3,m,kk;
int num_acc;
int atten;
int low1, mid1, high1, low2, mid2, high2, low3, mid3, high3;
double lon, lat, mw;
double rec_p[4];
double sigma;
double t_lon,t_lat,a10,b;
double a_est1, a_est2, a_est3;
double* a;
double c[5];
double* prob;
double* pro;
double result;
clock_t start,end;
rec_p[1] = 1./return_p[1]; // Return Period
rec_p[2] = 1./return_p[2]; // Return Period
rec_p[3] = 1./return_p[3]; // Return Period
///////////////////////////////////////////////////////////////////////
// Part 1 : Produce Seismic Hazard Map ///////////////////
///////////////////////////////////////////////////////////////////////
atten = 1; // Select Attenuation Relationship
///////////////////////////////////////////////////////////////////////
////////////////////////////////// Determine Attenuation Relationship /
if(atten == 1)
{
c[1] = 0.4854; c[2] = 1.2; c[3] = -0.8416; c[4] = -0.0061;
sigma = 0.8036;
}
else if(atten == 2)
{
c[1] = 0.5577; c[2] = 1.2; c[3] = -0.8587; c[4] = -0.0062;
sigma = 0.7629;
}
else if(atten == 3)
{
c[1] = 5.0244; c[2] = 0.5442; c[3] = -1.0020; c[4] = 0.0;
sigma = 0.1;
}
///////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////
////////////////////////////////// Set Target Acceleration ////////////
a = new double[33];
a[1]=0.001;a[2]=0.01;a[3]=0.02;a[4]=0.03;a[5]=0.04;a[6]=0.05;a[7]=0.06;a[8]=0.07;a[9]=0.08;a[10]=0.09;
a[11]=0.10;a[12]=0.11;a[13]=0.12;a[14]=0.13;a[15]=0.14;a[16]=0.15;a[17]=0.16;a[18]=0.17;a[19]=0.18;a[20]=0.19;
a[21]=0.20;a[22]=0.21;a[23]=0.22;a[24]=0.23;a[25]=0.24;a[26]=0.25;a[27]=0.26;a[28]=0.27;a[29]=0.28;a[30]=0.29;
a[31]=0.30;a[32]=0.31;
num_acc = 32;
///////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////
/////////////////////////////////// Read Gutenberg-Richter Param. a, b /
FILE* fin;
FILE* fout1;
FILE* fout2;
FILE* fout3;
FILE* fpro;
fin = fopen("ab_data_ceus.dat","r");
fout1 = fopen("seismic_map_500.dat","w");
fout2 = fopen("seismic_map_1000.dat","w");
fout3 = fopen("seismic_map_2400.dat","w");
for(j=1;j<=num_lat;j++)
{
for(i=1;i<=num_lon;i++)
{
m = (j-1)*num_lon + i;
fscanf(fin," %lf %lf %lf %lf %lf\n",&t_lon,&t_lat,&a10,&b,&mw);
aa[i][j] = a10;
bb[i][j] = b;
}
}
fclose(fin);
////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////
/////////////////////////////////// Read Normaly Distributed Function //
pro = new double[400];
fpro = fopen("CDF.txt","r");
for(i=0;i<=399;i++)
{
fscanf(fpro,"%lf",&pro[i]);
}
fclose(fpro);
////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////// Determine Acceleration of Target Return Period //////
prob = new double[17];
lat = S;
for(j=1;j<=num_lat;j++)
{
start = clock(); // check time
lon = W;
for(i=1;i<=num_lon;i++)
{
m = (j-1)*num_lon + i;
a_est1 = 0.;
a_est2 = 0.;
a_est3 = 0.;
kk = int(log(num_acc)/log(2));
low1 = low2 = low3 = 0;
high1 = high2 = high3 = num_acc;
mid1 = mid2 = mid3 = num_acc/2;
for(ii=1;ii<=num_acc;ii++)
prob[ii] = -1;
for(ii=1;ii<=kk;ii++)
{
k1 = mid1;
if( prob[k1] < 0 )
prob[k1] = probability(c,a[k1],sigma,lat,lon,pro);
if( prob[k1] < rec_p[1] ) { high1 = mid1; mid1 = (low1 + high1)/2; }
else if( prob[k1] > rec_p[1] ) { low1 = mid1; mid1 = (low1 + high1)/2; }
else { a_est1 = a[k1]; break; }
k2 = mid2;
if( prob[k2] < 0 )
prob[k2] = probability(c,a[k2],sigma,lat,lon,pro);
if( prob[k2] < rec_p[2] ) { high2 = mid2; mid2 = (low2 + high2)/2; }
else if( prob[k2] > rec_p[2] ) { low2 = mid2; mid2 = (low2 + high2)/2; }
else { a_est2 = a[k2]; break; }
k3 = mid3;
if( prob[k3] < 0 )
prob[k3] = probability(c,a[k3],sigma,lat,lon,pro);
if( prob[k3] < rec_p[3] ) { high3 = mid3; mid3 = (low3 + high3)/2; }
else if( prob[k3] > rec_p[3] ) { low3 = mid3; mid3 = (low3 + high3)/2; }
else { a_est3 = a[k3]; break; }
}
if(rec_p[1] > prob[k1])
{ if(k1 == 1) a_est1 = 0; a_est1 = a[k1] - (a[k1]-a[k1-1])*(prob[k1]-rec_p[1])/(prob[k1]-prob[k1-1]); }
else if(rec_p[1] < prob[k1])
a_est1 = a[k1] - (a[k1+1]-a[k1])*(prob[k1]-rec_p[1])/(prob[k1+1]-prob[k1]);
fprintf(fout1,"%.8f ",a_est1);
if(rec_p[2] > prob[k2])
{ if(k2 == 1) a_est2 = 0; a_est2 = a[k2] - (a[k2]-a[k2-1])*(prob[k2]-rec_p[2])/(prob[k2]-prob[k2-1]); }
else if(rec_p[2] < prob[k2])
a_est2 = a[k2] - (a[k2+1]-a[k2])*(prob[k2]-rec_p[2])/(prob[k2+1]-prob[k2]);
fprintf(fout2,"%.8f ",a_est2);
if(rec_p[3] > prob[k3])
{ if(k3 == 1) a_est3 = 0; a_est3 = a[k3] - (a[k3]-a[k3-1])*(prob[k3]-rec_p[3])/(prob[k3]-prob[k3-1]); }
else if(rec_p[3] < prob[k3])
a_est3 = a[k3] - (a[k3+1]-a[k3])*(prob[k3]-rec_p[3])/(prob[k3+1]-prob[k3]);
fprintf(fout3,"%.8f ",a_est3);
lon = lon + dD;
}
end = clock(); // check time
result=(double)(end-start) / CLOCKS_PER_SEC;
printf("%3.2f %2.1f sec %d\n",lat,result,iter);
fprintf(fout1,"\n");
fprintf(fout2,"\n");
fprintf(fout3,"\n");
lat = lat + dD;
}
fclose(fout1);
fclose(fout2);
fclose(fout3);
//////////////////////////////////////////////////////////////////////////
return 0;
}
//////////////////////////////////////////////////////////////////////////////
// Yearly generate probability above u0 acceleration //
//////////////////////////////////////////////////////////////////////////////
double probability(double c[5],double u,double sigma,double lat2,double lon2,double* pro)
{
int i, j, k, l;
int f_M_check;
double lat, lon, M;
double alpha, beta;
double pi, dv, dh, z;
double CDF, p_M, f_MM;
double a, R;
double prob;
double* f_M;
f_M = new double[num_M+1];
pi = 3.141592;
dv = 2*pi*6400/360; // unit distance of latitude
prob = 0;
f_M_check = 1;
lon = W;
for(i=1;i<=num_lon;i++)
{
lat = S;
for(j=1;j<=num_lat;j++)
{
beta = 2.303*bb[i][j];
if(aa[i][j] != 0)
{
alpha = 2.303*log10(aa[i][j]);
p_M = exp(alpha-beta*Ml); // V
dh = dv*cos((lat+lat2+dD)/2*pi/180); // unit distance of longitude
R = sqrt(pow(((lon2-lon)*dh),2)+pow(((lat2-lat)*dv),2));
M = Ml;
for(k=1;k<=num_M;k++)
{
if( control_b == 10 )
{
f_MM = beta*exp(-beta*(M+dM/2-Ml))/(1-exp(-beta*(Mu-Ml)))*dM; // P[M]
f_M[k] = f_MM;
}
else
{
if( f_M_check != 0 )
{
f_MM = beta*exp(-beta*(M+dM/2-Ml))/(1-exp(-beta*(Mu-Ml)))*dM; // P[M]
f_M[k] = f_MM;
if( k == num_M ) f_M_check = 0;
}
}
a = exp(c[1]+c[2]*(M+dM/2)+c[3]*log(sqrt(R*R+100))+c[4]*(sqrt(R*R+100)))/981; // Acc. by ARship
if( a > 0.001 )
{
z = (log(u*981)-log(a*981))/sigma;
l = int(ceil(fabs(z)*100-0.5));
if(l > 399) l=399;
if(z < 0)
CDF = 1-pro[l];
else
CDF = pro[l];
prob = prob + f_M[k]*CDF*p_M; // Annual Occurence Rate
}
M=M+dM;
}
}
lat = lat + dD;
}
lon = lon + dD;
}
return prob;
}
