🚀 go-pugleaf

zelfs als ik dat aanpas gaat het nog niet;...ik zal es de rest van mijn
programma meesturen, maar voorlopig is het enkel het inlezen van de data dat
belangrijk is...zolang hij dat niet doet kan ik toch ni verder....

--
------------------------------------
Matthias Demuzere
O.L. Vrouwstraat 50
3000 Leuven

<Bart.Vanhauwaert@nowhere.be> schreef in bericht
news:u3qmta.3e1.ln@ns1.irule.be...
> Matthias Demuzere <Matthias.demuzere@student.kuleuven.ac.be> wrote:
> > Ik heb het volgende probleem. Als ik dit stukje programma laat lopen,
dan
> > geeft hij een waarde aan mn j van 4067520 (of toch zoiets in die grootte
> > orde...) Iemand aan idee waar hij die waarde vandaan haalt???
> > alvast bedankt
>
> > Matthias
>
> > -------------------------------------------
> > int j;
>
> > FILE *ifp;
> >   FILE *ofp;
>
> >   /*reading file*/
>
> >   ifp = fopen("../gjelbreen","r");
>
> >   /*reading data from file*/
>
> >   for (j=1;j<=MAXX;j++)
>
> >   {
>
> >    fscanf(ifp,"%lf%lf%lf",&x[j],&hb[j],&s[j]);
> >    printf("%f",hb[j]);
>
> >   }
>
> Zonder het volledige programma te zien, vind ik je for loop toch een
> beetje twijfelachtig. Normaal telt men in C van 0 tot MAXX-1 en niet van
> 1 tot MAXX. Dus schrijf je
>
> for (j=0; j<MAXX; ++j)
> ...
>
>
> Als je de variable x bv als
>
> double x[MAXX];
>
> hebt gedeclareerd, dan ga je dus een waarde schrijven buiten de limieten
> van je array. En het kan wel eens heel goed dat op die plaats j staat en
> je daarom dus onvoorziene toestanden krijgt.
>
> cu bart
>
> --
> http://www.irule.be/bvh/

zelfs als ik dat aanpas gaat het nog niet;...ik zal es de rest van mijn
programma meesturen, maar voorlopig is het enkel het inlezen van de data dat
belangrijk is...zolang hij dat niet doet kan ik toch ni verder....

/*standard libraries*/

#include <stdio.h>  /*standard input-output*/
#include <stdlib.h> /*standard C library*/
#include <math.h>  /*math library*/


/*definieren van de parameters*/


#define GA   9.81      /*gravity acceleration (m/s�)*/
#define K       (2.1*SEC_PER_YEAR)       /*heat conductivity (W/mK)*/
#define RHO      910.0        /*ice density (kg/m�)*/
#define G       (0.042*SEC_PER_YEAR)     /*geothermal heat flux (*/
#define CP       2009.0       /*specific heat capacity (J/kgK)*/
#define TABS     273.15       /* absolute temperature (K)*/
#define TSEA                  /*temperature at sea level (K)*/
#define BETA  0.000000098     /*constant parameter (K/Pa)*/
#define ADIAT     (-0.012)       /*Factor regarding temperature/height
(K/m)*/
#define A        1.0e-16                     /*flow parameter */

#define TIMETOT  50000.0       /*total calculation time (year)*/
#define MAXX  122       /*maximum number of elements*/
#define MAXXARRAY   (MAXX+1)
#define MAXZ  51
#define MAXZARRAY (MAXZ+1)
#define SEC_PER_YEAR 31556926     /*numbers of seconds in 1 year*/
#define DT   10.0
#define DZETA  0.02      /*interval on the vertical axes*/
#define DX    10.000      /*interval on the horizontal axes (m)*/




/*main program*/

int main(int argc, char **argv)   /*arguments to communicate with the main
system*/

 {
  int i,j,k,tlapse;
  double t0[MAXZARRAY][MAXXARRAY], t1[MAXZARRAY][MAXXARRAY],
w[MAXZARRAY][MAXXARRAY], u[MAXZARRAY][MAXXARRAY];
  double h[MAXXARRAY];      /*ice sheet thickness*/
  double s[MAXXARRAY];      /*Height of surface*/
  double hb[MAXXARRAY];      /*height of basis of glacier*/
  double x[MAXXARRAY];
  double m[MAXXARRAY];
  double zeta[MAXZARRAY];
  double tsurf[MAXXARRAY];
  double tpmp[MAXZARRAY][MAXXARRAY];
  double wsurf[MAXXARRAY];
  double f[MAXZARRAY];
  double g[MAXZARRAY];
  double alfa[MAXZARRAY];
  double beta[MAXZARRAY];
  double gamma[MAXZARRAY];
  double delta[MAXZARRAY];
  double taud[MAXXARRAY];

  FILE *ifp;
  FILE *ofp;

  /*reading file*/

  ifp = fopen("../gjelbreen","r");

  /*reading data from file*/

  for (j=0;j<MAXX-1;j++)

  {

   fscanf(ifp,"%lf%lf%lf",&x[j],&hb[j],&s[j]);
   printf("%f",hb[j]);

  }


exit(1);

----------------------------------------------------------------------------
-------------------------------------

/*constructing matrix,going from top to bottom, evolving downstream*/

  /*
  initializing temperature and velocity profile
  i=1: top of de ice column at height H
  i=MAXZ: bottom of the ice column at height 0
  j=1: most upper part of the ice mass on the horizontal axes
  j=MAXX: lowest part of the ice mass on the horizontal axes
  */


  for(j=2;j<=MAXX-1;j++)

  {

   /*parameters evolving in x-direction*/

   h[j]=s[j]-hb[j];

   /*calculating mass balance and temperature at sea surface */

   if (s[j]<1000)

   {

    m[j]=0.50-(double)(0.00025*hb[j]);  /*mass balance is lineair with
elevation < 1000m*/
    tsurf[j]=257.4-(double)(0.0051*hb[j]); /*TSURF is lineair with elevation
< 1000m*/

   }

   else

   {
    m[j]=(double)(0.678*exp((-0.00106)*hb[j]));   /*exponential decrease of
the mass balance for elevation > 1000m*/
    tsurf[j]=266.4-(double)(0.0143*hb[j]); /*TSURF is lineair with elevation
> 1000m*/

   }

   /*vertical speed at surface*/

   wsurf[j]=(double)(-m[j]);

   /*driving stress used in equation for horizontal speed*/

   taud[j]=(double)(-RHO*GA*h[j]*((s[j+1]-s[j-1])/2*DX));



   for(i=1;i<=MAXZ;i++)

   {

    /*parameters evolving in "z"-direction*/

    zeta[i]=(double)((i-1)*DZETA);

    w[i][j]=wsurf[j]*(1+(double)((i-1)*DZETA));
    u[i][j]=((1.0-pow(zeta[i],4.0))*((2.0*A*h[j])/5.0)*pow(taud[j],3));


   }

   for(i=1;i<=MAXZ;i++)

   {

    /*boundary conditions at ice divide*/

    t1[i][1]=233.15;
    u[i][1]=-u[i][3]; /*at point i,1, u is measured in the opposite
direction*/
    w[i][1]=w[i][3]; /*the vertical speed is in the same direction*/
   }

  }


  /*

  time loop

  */

  tlapse =(int)(TIMETOT/DT+0.5);

  for(k=1;k<=tlapse;k++)

  {

   /*temperature calculation*/

   for(j=2;j<=MAXX-1;j++)

   {

    /*define f,g at top*/

    f[1]=0;
    g[1]=TSEA+(ADIAT*s[j]);

    for(i=2;i<MAXZ-1;i++)

    {
     /*calculate alfa,beta,gamma,delta to solve the implicit equation*/

     alfa[i]=-(K*DT/RHO*CP*pow(h[j],2)*pow(DZETA,2))+(w[i][j]/
       h[j]*DZETA)-(u[i][j]*(s[j+1]-s[j-1])/(2.0*h[j]*DX*DZETA))-
       (u[i][j]*zeta[i]*(h[j+1]-h[j-1])/(2.0*DX*DZETA*h[j]));

     beta[i]=1.0+(K*DT/(RHO*CP*pow(h[j],2)*pow(DZETA,2)))-
       (w[i][j]/h[j]*DZETA)+(u[i][j]*(s[j+1]-s[j-1])/2.0*h[j]*DX*DZETA)
       +(u[i][j]*zeta[i]*(h[j+1]-h[j-1])/2.0*h[j]*DZETA*DX)+(u[i][j]/DX);

     gamma[i]=K*DT/CP*pow(h[j],2)*pow(DZETA,2);

     delta[i]=t0[i][j]+(zeta[i]*GA*(s[j+1]-s[j-1])*(u[i][j]-u[i-1][j])/
        CP*2.0*DX*DZETA)+((u[i][j]/DX)*t0[i][j-1]);


     f[i]=((-gamma[i])/((alfa[i]*f[i-1])+beta[i]));
     g[i]=(delta[i]-(alfa[i]*g[i-1])/(alfa[i]*f[i-1]+beta[i]));

     /*calculating temperature*/

     t0[i][j]=f[i]*t1[i][j]+g[i];


     /*explicit method*/


/*t1[i][j]=t0[i][j]+((K*DT)/(RHO*CP*H[j]*H[j]))((t0[i-1][j-1]-(2*t0[i][j])+t
0[i+1][j+1])

/(DZETA*DZETA))+((t0[i][j]-t0[i-1][j])/DZETA)((w[i][j]/H[j])-(u[i][j]*

((1.0/H[j])*((S[j+1]-S[j-1])/(2.0*DX)-(zeta[i]*((H[j+1]-H[j-1])/(2.0*DX)))))
)
         -u[i][j]*((t0[i][j]-t1[i][j-1])/DX)+(zeta[i]*GA*(S[j+1]-S[j-1])*(u[
i][j]-u[i-1][j]))
         /(CP*2.0*DX*DZETA));*/

     /*T can't be higher than pressure melting point*/

     tpmp[i][j]=TABS-(double)(BETA*RHO*GA*h[j]*zeta[i]);

     if (t1[i][j]>tpmp[i][j])

     t1[i][j]=tpmp[i][j];

    }

    /*calculate boundary conditions for t and velocity at the end of the ice
mass*/

    for (i=MAXZ-1;i>=1;i--)

    {
     t1[i][MAXX]=t1[i][MAXX-1];
     u[i][MAXX]=u[i][MAXX-1];
     w[i][MAXX]=w[i][MAXX-1];

    }


    /*boundary conditions*/

    t1[1][j]=tsurf[j];
    t1[MAXZ][j]=(DZETA*G*h[j]/K)+t0[MAXZ-1][j];

   }

   /*at every time step,basal temperatures must be displayed on the screen*/

   printf("Time=%.0f year --- Basal temp=%fC\n",((float)(k-1)*DT),
      (t1[MAXZ][MAXX]-TABS));

   /*replacing old temperature at t by new at t+1*/

   for (j=1;j<=MAXX;j++)

   {
    for (i=1;i<=MAXZ;i++)

    {

     t0[i][j]=t1[i][j];

    }

   }

  }

  /*writing the output on file outfp*/

  ofp=fopen("tprofout.xls","w");

  for (j=1;j<=MAXX;j++)

  {

   for (i=1;i<=MAXZ;i++)

   {

    fprintf(ofp,"%f\t%f\t%f\n%f\n",x[j],t0[i][j],w[i][j],u[i][j]);
   }

  }

  fclose(ofp);
  return (0);

 }