Subject: keywords and parameters for the farside pipeline From: Irene Gonzalez-Hernandez Date: Tue, 28 Jun 2005 14:29:14 -0700 To: pscherrer@solar.stanford.edu CC: IRENEGH@noao.edu, John Leibacher Dear Phil and Kenny, We keep two independent programs to do the postel projection of the input images, called map_gong.c and map_mdi.c. I'm attaching the version of map_mdi.c in case it maybe helpful to you. But here are the parameters and keyword values that map_mdi uses: Parameters: For the rotation correction and truncation: #define ROT_STATIST_MARGIN 0.100 /* Solar radii (inside of limb) */ #define TRUNC_MARGIN 0.004 /* Solar radii (inside of limb) */ #define TRUNC_MARGIN_PIX 1.5 /* pixel (inside of limb) */ For the definition of the origin of coordinates: #define XCORRECT 1.0 /* For MDI center of lower left */ #define YCORRECT 1.0 /* pix is (1, 1) --> (0, 0) */ keyword values: FNDLMBXC FNDLMBYC R_SUN FNDLMBAN OBS_L0 OBS_B0 SOLAR_P OBS_R0 (Default 1 AU) PLATE_X PLATE_Y DPC_CROP (Defaulted to (1. - TRUNC_MARGIN)*sun_rad_min) ** This keyword exists in the full res (1024x1024) MDI images but not in low res. T_STEP (Default 60s) T_OBS (Default 1996.01.01_00:00:00_TAI) ** I'm not quite sure why there is a default value for this one ?? For the low resolution (192x192) MDI images, we also use a preprocesor that separates the cubes of images into individual ones. It reads from the cube header the keywords below for each slide and put them into the individual image headers. {"T_OBS", "%s"}, {"OBS_L0", "%e"}, {"OBS_B0", "%e"}, {"OBS_R0", "%e"}, {"SOLAR_P", "%e"}, {"FNDLMBXC", "%e"}, {"FNDLMBYC", "%e"}, {"R_SUN", "%e"}, {"FNDLMBAN", "%e"}, {"PLATE_X", "%e"}, {"PLATE_Y", "%e"}, {"SOLAR_P", "%e"}, {"QUALITY", "%x"}, {"DATAMIN", "%e"}, {"DATAMAX", "%e"} Hope this helps, cheers, Irene. /* Program map_mdi * This program processes a SOHO MDI image for helioseismic computations. * * It performs the functions of three separate algorithms: * * (1) map_perspect, which maps the image onto a representation equivalent * the perspective from a position not necessarily the same as that * viewed by the SOHO spacecraft, * * (2) derotate_mdi_image, which measures and removes solar rotation * before the mapping (1) is done, and * * (3) Padding of the margins of the image, greatly reducing the ringing * of the Fourier interpolator in (1) caused by the discontinuity * across the solar limb. * * The original version of this program was written by Seth Redfield * in the summer of 1995. It was revised by Mark Fagan in the summer of * 1996 to incorporate the option of inputing a list of perspectives and * corresponding output filenames. It has the additional flexibility of * mapping the center of the projection anywhere on the map, as specified * by program arguments 6 and 7. The optional "mode" parameter specifies * the geometry of the projection as follows: * * -los (default) Line-of-sight projection onto the data * array through an imaginary pinhole. * * -perp Perpendicular projection of the Sun's * spherical surface onto the data array. * * -postel Postel projection. * * Compile this program with the command * * "gcc -O3 map_mdi.c -o map_mdi -lm". */ #include #include #include "/home/irenegh/src/fspipe/fspipe_2.0/sources/do_fits.c" #define LONG_LINE_LENGTH 160 #define WRONG_NUM_COLS -1 #define ROW_FAILURE -2 #define FORWARD 1 #define INVERSE -1 #define ROT_STATIST_MARGIN 0.100 /* Solar radii (inside of limb) */ #define TRUNC_MARGIN 0.004 /* Solar radii (inside of limb) */ #define TRUNC_MARGIN_PIX 1.5 /* pixel (inside of limb) */ #define XCORRECT 1.0 /* For MDI center of lower left */ #define YCORRECT 1.0 /* pix is (1, 1) --> (0, 0) */ /* #define LIMB_TRUNC_VALUE -400.0 */ /* Diagnostic */ /* #define DARK_SIDE_VALUE -1000.0 */ /* Diagnostic */ #define SPACE_VALUE 0.0 #define LIMB_TRUNC_VALUE 0.0 #define DARK_SIDE_VALUE 0.0 #define DEFAULT_PARITY 1 #define LOS 0 #define PERP 1 #define POSTEL 2 #define PI 3.14159265359 #define DEFAULT_P_ANGLE -0.17 /* degree */ #define DEFAULT_SUN_RADIUS asin(695997./149597870.) /* radian */ int print_diag(int argc); int get_params(FILE *fp, char *argv[], float *plon_ref, float *plat_ref, float *px0_out, float *py0_out, char out_file[]); int lp2ge(int arg); int inset_data(float buff[], int nx1, int nx2, int ny1, int ny2); int expand_frame(float buff[]); int read_data(int fid, int n_dat, int n_buff, float **pbuff); int interp_data(float *buff, int magarg); int f_interp(int num_pts, float dat[]); int rfft(int n, float dat_in[], float dat_out[], int mode); int count_bits(int npts); int fft(int npts, float dat[], int mode); int bit_rev(int npts, float dat[]); int perp_project(float lon_ref, float lat_ref, float buff_in[], float buff_out[]); int los_project(float lon_ref, float lat_ref, float buff_in[], float buff_out[]); int postel(float lon_ref, float lat_ref, float buff_in[], float buff_out[]); float bicubic_interp(float x, float y, float **dat_in); float cu_interp(float ref_grid[], float x); int gile(float x); int pad_image(float data[], float xc, float yc, float rad); int truncate_limb(float data_in[], float xc, float yc, float rad, float value); int meas_rot(float data[], float *prefx, float *prefy, float *pmean, float *pgradx, float *pgrady); float remove_rot(float data[], float refx, float refy, float mean, float gradx, float grady); int getdt(char line[], float *pdtime); int gettime(char line[], int *ptime); int date_to_jday(int year, int month, int dom); int jday_to_date(int day, char alpha_date[], int alpha_set); int isleap(int year); static float pi = PI; static float rad_fact = PI/180.; static int nxin, nyin, nxout, nyout; static int nxframe, nyframe, mag, nxmagframe, nymagframe; static int parity; static float dxmagframe, dymagframe, dxout, dyout; static float x_center, y_center; static float x_center_mag, y_center_mag, x_center_out, y_center_out; static float dtime, b0, l0, pangle; static float sun_rad_radian, sun_rad_arcsec; static float sun_rad_maj, sun_rad_min, lon_ref, lat_ref; static float maj_ax_fact, min_ax_fact, maj_ax_ang; main(int argc, char *argv[]) { char line[LINE_LENGTH + 1]; char alpha_date[LINE_LENGTH + 1]; char alpha_time[LINE_LENGTH + 1]; char out_file[LINE_LENGTH + 1]; int fid, fd, stat, mode; int dt_stat; int numbits, diag, num_params, *params; int time_secs, soho_day, day_secs, hrs, mins, secs; float rms_dev; float *data_in, *data_out; float refx, refy, mean, gradx, grady; float trunc_rad_pix, trunc_xc, trunc_yc; FILE *fp, *fopen(); init_fits(argv[0]); if (argc != 4 && argc != 5 && argc != 8 && argc != 9) print_diag(argc); mode = LOS; if (argc == 5 || argc == 9) { if (seg_cmp(argv[argc - 1], "-los", 5) == 0) mode = LOS; else if (seg_cmp(argv[argc - 1], "-perp", 6) == 0) mode = PERP; else if (seg_cmp(argv[argc - 1], "-postel", 8) == 0) mode = POSTEL; else print_diag(argc); } parity = DEFAULT_PARITY; mag = 2; /* Open the reference file and read the fits head. */ if ((fid = open(argv[2], O_RDONLY)) == EOF) { fprintf(stderr, "%s: Cannot open file named %s. Exiting.\n", dfprogname, argv[2]); exit(0); } if (get_head(fid, &inhead) == EOF) { fprintf(stderr, "%s: Defective fits head in file %s. Exiting.\n", dfprogname, argv[2]); exit(0); } if (read_head(inhead) == EOF) { fprintf(stderr, "%s: Cannot read fits head in file %s. Exiting.\n", dfprogname, argv[2]); exit(0); } close(fid); nxout = nx; nyout = ny; dxout = dx; dyout = dy; data_out = (float *)malloc(nxout*nyout*sizeof(float)); if (naxis != 2) { fprintf(stderr, "%s: File %s: Array must be 2 dimensions. ", dfprogname, argv[1]); fprintf(stderr, "It is %d. Exiting.\n", naxis); exit(0); } /* Open the input file and read the fits head. */ if ((fid = open(argv[1], O_RDONLY)) == EOF) { fprintf(stderr, "%s: Cannot open file named %s. Exiting.\n", dfprogname, argv[1]); exit(1); } if (get_head(fid, &inhead) == EOF) { fprintf(stderr, "%s: Defective fits head in file %s. Exiting.\n", dfprogname, argv[1]); exit(1); } if (read_head(inhead) == EOF) { fprintf(stderr, "%s: Cannot read fits head in file %s. Exiting.\n", dfprogname, argv[1]); exit(1); } if (naxis != 2) { fprintf(stderr, "%s: File %s: Array must be 2 dimensions. ", dfprogname, argv[1]); fprintf(stderr, "It is %d. Exiting.\n", naxis); exit(1); } /* Read needed keyword values from FITS head. */ if (get_keyword_value(inhead, "FNDLMBXC", "%e", &x_center) == EOF) { fprintf(stderr, "%s: Keyword \"FNDLMBXC\" missing in file %s. ", dfprogname, argv[1]); fprintf(stderr, "Exiting.\n"); exit(1); } x_center -= XCORRECT; if (get_keyword_value(inhead, "FNDLMBYC", "%e", &y_center) == EOF) { fprintf(stderr, "%s: Keyword \"FNDLMBYC\" missing in file %s. ", dfprogname, argv[1]); fprintf(stderr, "Exiting.\n"); exit(1); } y_center -= YCORRECT; if (get_keyword_value(inhead, "R_SUN", "%e", &sun_rad_min) == EOF) { fprintf(stderr, "%s: Keyword \"R_SUN\" missing in file %s. ", dfprogname, argv[1]); fprintf(stderr, "Exiting.\n"); exit(1); } sun_rad_maj = sun_rad_min; if (get_keyword_value(inhead, "FNDLMBAN", "%e", &maj_ax_ang) == EOF) { fprintf(stderr, "%s: Keyword \"FNDLMBAN\" missing in file %s. ", dfprogname, argv[1]); fprintf(stderr, "Exiting.\n"); exit(1); } if (get_keyword_value(inhead, "OBS_L0", "%e", &l0) == EOF) { fprintf(stderr, "%s: Keyword \"OBS_L0\" missing in file %s. ", dfprogname, argv[1]); fprintf(stderr, "Exiting.\n"); exit(1); } if (get_keyword_value(inhead, "OBS_B0", "%e", &b0) == EOF) { fprintf(stderr, "%s: Keyword \"OBS_B0\" missing in file %s. ", dfprogname, argv[1]); fprintf(stderr, "Exiting.\n"); exit(1); } /* if (get_keyword_value(inhead, "SOLAR_P", "%e", &pangle) == EOF) { fprintf(stderr, "%s: Keyword \"SOLAR_P\" missing in file %s. ", dfprogname, argv[1]); fprintf(stderr, "Defaulting to %.1f.\n", DEFAULT_P_ANGLE); pangle = DEFAULT_P_ANGLE; } */ /* Commented out because in lev1.8 SOLAR_P = inaccurate 0.0 */ pangle = DEFAULT_P_ANGLE; if (get_keyword_value(inhead, "OBS_R0", "%e", &sun_rad_arcsec) == EOF) { fprintf(stderr, "%s: Keyword \"OBS_R0\" missing in file %s.\n", dfprogname, argv[1]); fprintf(stderr, "Defaulting to angular radius of Sun at 1 AU "); fprintf(stderr, "= %e arcsec.\n", DEFAULT_SUN_RADIUS); sun_rad_arcsec = DEFAULT_SUN_RADIUS; sun_rad_radian = sun_rad_arcsec*pi/(3600.*180.); } sun_rad_radian = DEFAULT_SUN_RADIUS; if (get_keyword_value(inhead, "PLATE_X", "%e", &trunc_xc) == EOF) { fprintf(stderr, "%s: Keyword \"PLATE_X\" missing in file %s.\n"); fprintf(stderr, "Defaulting to center: .5*(NAXIS1 - 1.) = %.1f\n", .5*(nx - 1.)); trunc_xc = .5*(nx - 1.); } if (get_keyword_value(inhead, "PLATE_Y", "%e", &trunc_yc) == EOF) { fprintf(stderr, "%s: Keyword \"PLATE_Y\" missing in file %s.\n"); fprintf(stderr, "Defaulting to center: .5*(NAXIS2 - 1.) = %.1f\n", .5*(ny - 1.)); trunc_xc = .5*(ny - 1.); } if (get_keyword_value(inhead, "DPC_CROP", "%e", &trunc_rad_pix) == EOF) { fprintf(stderr, "%s: Keyword \"DPC_CROP\" missing in file %s. "); fprintf(stderr, "Defaulting to %f solar radius.\n", TRUNC_MARGIN); trunc_rad_pix = (1. - TRUNC_MARGIN)*sun_rad_min; } trunc_rad_pix -= TRUNC_MARGIN_PIX; if ((dt_stat = get_keyword_value(inhead, "T_STEP", "%s", line)) == EOF) { fprintf(stderr, "%s: Keyword T_STEP missing in file %s.\n", dfprogname, argv[1]); fprintf(stderr, "Defaulting to 60 sec.\n"); sprintf(line, "'60s'"); } getdt(line, &dtime); if (get_keyword_value(inhead, "T_OBS", "%s", line) == EOF) { fprintf(stderr, "%s: Keyword \"T_OBS\" missing in file %s. ", dfprogname, argv[1]); fprintf(stderr, "Defaulting to 1996.0.\n"); sprintf(line, "'1996.01.01_00:00:00_TAI'"); } /* Use the alphanumeric time to compute the number of seconds since * 1996.0. */ gettime(line, &time_secs); /* Set dimensions and spacings for magnified frame. */ nxin = nx; nyin = ny; nxframe = lp2ge(nx); nyframe = lp2ge(ny); nxmagframe = mag*nxframe; nymagframe = mag*nyframe; dxmagframe = dx/(float)mag; dymagframe = dy/(float)mag; /* Now read the data in the input file into a magnified array. */ read_data(fid, nxin*nyin, nxmagframe*nymagframe, &data_in); close (fid); if (bitpix < 0) bitpix = -bitpix; /* Optional diagnostic. */ /* fprintf(stdin, "\n%s: WRITING DIAGNOSTIC FILE DIAG1.\n"); fid = open("DIAG1", O_CREAT | O_WRONLY | O_TRUNC, 0644); current_head = NULL; write_fits(fid, data_in); close(fid); */ /* Measure the mean Doppler shift and rotation gradient and remove * them from the image. For flexibility, the following two functions * should (and do) have programs to themselves. They are included * here for efficiency, to save the need to read and write the file * one additional time for each. */ if (meas_rot(data_in, &refx, &refy, &mean, &gradx, &grady) == EOF) { fprintf(stderr, "%s: Cannot measure rotation gradient. Exiting.\n", dfprogname); exit(0); } rms_dev = remove_rot(data_in, refx, refy, mean, gradx, grady); fprintf(stderr, "\n%s: Diagnostics for file %s:\n", dfprogname, argv[1]); fprintf(stderr, "%11.4e %11.4e %11.4e %11.4e %11.4e %11.4e\n", refx, refy, mean, gradx, grady, rms_dev); /* Optional diagnostic. */ /* fprintf(stdin, "\n%s: WRITING DIAGNOSTIC FILE DIAG2.\n"); fid = open("DIAG2", O_CREAT | O_WRONLY | O_TRUNC, 0644); current_head = NULL; write_fits(fid, data_in); close(fid); */ /* Expand the data to occupy a larger 2-dimensional array, (nxframe, nyframe) if this is larger than (nxin, nyin) */ inset_data(data_in, nxin, nxframe, nyin, nyframe); nx = nxframe; ny = nyframe; /* Pad the margins of the image so that the Fourier magnifier will * not ring on the limb discontinuities. */ pad_image(data_in, trunc_xc, trunc_yc, trunc_rad_pix); /* Optional diagnostic. */ /* fprintf(stdin, "\n%s: WRITING DIAGNOSTIC FILE DIAG3.\n"); fid = open("DIAG3", O_CREAT | O_WRONLY | O_TRUNC, 0644); current_head = NULL; write_fits(fid, data_in); close(fid); */ /* Expand the data to occupy the larger 2-dimensional array, * (nxmagframe, nymagframe), which is larger than * (nxframe, nyframe). */ inset_data(data_in, nxframe, nxmagframe, nyframe, nymagframe); nx = nxmagframe; ny = nymagframe; /* Redistribute the data points onto the expanded scale. */ expand_frame(data_in); sun_rad_maj *= (float)mag; sun_rad_min *= (float)mag; x_center_mag = x_center*mag; y_center_mag = y_center*mag; nx = nxmagframe; ny = nymagframe; dx = dxmagframe; dy = dymagframe; /* Interpolate the data between the data points in the expanded array. */ interp_data(data_in, mag); /* Optional diagnostic. */ /* fprintf(stdin, "\n%s: WRITING DIAGNOSTIC FILE DIAG4.\n"); fid = open("DIAG4", O_CREAT | O_WRONLY | O_TRUNC, 0644); current_head = NULL; write_fits(fid, data_in); close(fid); */ truncate_limb(data_in, mag*trunc_xc, mag*trunc_yc, mag*trunc_rad_pix, LIMB_TRUNC_VALUE); /* Optional diagnostic. */ /* fprintf(stdin, "\n%s: WRITING DIAGNOSTIC FILE DIAG5.\n"); fid = open("DIAG5", O_CREAT | O_WRONLY | O_TRUNC, 0644); current_head = NULL; write_fits(fid, data_in); close(fid); */ /* Open list file and loop through lon_ref, lat_ref, and out_file */ fp = NULL; if ((argc == 4 || argc == 5) && (fp = fopen(argv[3], "r")) == NULL) { fprintf(stderr, "%s: Cannot open file named %s. Exiting.\n", dfprogname, argv[3]); exit(0); } printf("\n"); while (get_params(fp, argv, &lon_ref, &lat_ref, &x_center_out, &y_center_out, out_file) != EOF) { fprintf(stderr, "%s: Writing file %s\n", dfprogname, out_file); nx = nxmagframe; ny = nymagframe; dx = dxmagframe; dy = dymagframe; x_center = x_center_mag; y_center = y_center_mag; if (mode == LOS) los_project(lon_ref, lat_ref, data_in, data_out); else if (mode == PERP) perp_project(lon_ref, lat_ref, data_in, data_out); else if (mode == POSTEL) postel(lon_ref, lat_ref, data_in, data_out); /* Write out the data. */ if ((fid = open(out_file, O_CREAT | O_EXCL | O_WRONLY, 0644)) == EOF) { fprintf(stderr, "%s: File named %s exists. Overwrite? (y/n): ", dfprogname, out_file); fgets(line, LINE_LENGTH + 1, stdin); if (seg_cmp(line, "y", 1) != 0) exit(0); if ((fid = open(out_file, O_CREAT | O_TRUNC | O_WRONLY, 0644)) == EOF) { fprintf(stderr, "%s: Cannot open file %s for write. Exiting.\n", dfprogname, out_file); exit(0); } } nx = nxout; ny = nyout; dx = dxout; dy = dyout; x_center = x_center_out; y_center = y_center_out; num_head_recs = make_head(¤t_head); insert_keyword_value(current_head, 0, "OBS_L0", "%e", l0, "/ degree"); insert_keyword_value(current_head, 0, "OBS_B0", "%e", b0, "/ degree"); insert_keyword_value(current_head, 0, "REF_L0", "%e", lon_ref, "/ degree"); insert_keyword_value(current_head, 0, "REF_B0", "%e", lat_ref, "/ degree"); insert_keyword_value(current_head, 0, "X0", "%e", x_center, "/ pixel"); insert_keyword_value(current_head, 0, "Y0", "%e", y_center, "/ pixel"); insert_keyword_value(current_head, 0, "OBS_R0", "%e", sun_rad_radian, "/ radian"); insert_keyword_value(current_head, 0, "T_STEP", "%e", dtime, "/ sec"); soho_day = time_secs/86400; day_secs = time_secs%86400; hrs = day_secs/3600; mins = (day_secs/60)%60; secs = (day_secs%3600)%60; jday_to_date(soho_day, alpha_date, 0); sprintf(alpha_time, "'%s %02d:%02d'", alpha_date, hrs, mins); insert_keyword_value(current_head, 0, "MDITIME", "%i", time_secs, "/ sec since 1996.0"); insert_keyword_value(current_head, 0, "CALTIME", "%s", alpha_time, "/ yr mo dy hr:min"); write_fits(fid, data_out); close(fid); } if (fp != NULL) fclose(fp); free(data_in); free(data_out); } int print_diag(int argc) { FILE *dest; if (argc == 1) dest = stdout; else dest = stderr; fprintf(dest, "\nUsage: %s in_file ref_file lon_ref lat_ref ", dfprogname); fprintf(dest, "x_cen y_cen out_file [mode]\n"); fprintf(dest, "\n %s in_file ref_file list_file [mode]\n", dfprogname); fprintf(dest, "\nIn the second case, list_file is an ascii "); fprintf(dest, "list of the form\n\n"); fprintf(dest, " lon_ref_1 lat_ref_1 x_cen_1 y_cen_1 "); fprintf(dest, "out_file_1\n"); fprintf(dest, " lon_ref_2 lat_ref_2 x_cen_2 y_cen_2 "); fprintf(dest, "out_file_2\n"); fprintf(dest, " ..... ..... ..... ..... .....\n\n"); fprintf(dest, "The units of longitude and latitude are degrees.\n\n"); fprintf(dest, "The string \"obs\" can be input "); fprintf(dest, "for lon_ref and/or lat_ref.\n\n"); fprintf(dest, "The units of x_cen and y_cen are pixels.\n\n"); fprintf(dest, "The string \"nom\" can be input "); fprintf(dest, "for x_cen and/or y_cen to center\n"); fprintf(dest, "the projection in the output frame.\n\n"); fprintf(dest, "The optional mode can take any of the "); fprintf(dest, "following three values:\n\n"); fprintf(dest, " -los (default) Line-of-sight "); fprintf(dest, "projection onto the data array\n"); fprintf(dest, " through an "); fprintf(dest, "imaginary pinhole.\n\n"); fprintf(dest, " -perp Perpendicular "); fprintf(dest, "projection of the Sun's spherical\n"); fprintf(dest, " surface onto "); fprintf(dest, "the data array.\n\n"); fprintf(dest, " -postel Postel projection.\n\n"); if (argc == 1) exit(0); exit(1); } /* Function get_params * This function looks up parameters (reference longitude and latitude, * x- and y-coordinates for center of image and output filename) for * the output files. It determines whether the parameters are expressed * in the program arguments or are in a file named by a program * argument, in which case it reads the parameters therefrom. */ int get_params(FILE *fp, char *argv[], float *plon_ref, float *plat_ref, float *px0_out, float *py0_out, char out_file[]) { char line[LONG_LINE_LENGTH + 1]; char lon_string[LINE_LENGTH + 1]; char lat_string[LINE_LENGTH + 1]; char x0_string[LINE_LENGTH + 1]; char y0_string[LINE_LENGTH + 1]; int lon_stat, lat_stat, x0_stat, y0_stat; static int start = 1; if (fp != NULL) /* Parameters contained in file to which fp points */ { while(fgets(line, LONG_LINE_LENGTH + 1, fp) != NULL) { if (sscanf(line, "%s %s %s %s %s", lon_string, lat_string, x0_string, y0_string, out_file) == 5) { if (seg_cmp(lon_string, "obs", 3) == 0) { *plon_ref = l0; lon_stat = 1; } else lon_stat = sscanf(lon_string, "%f", plon_ref); if (seg_cmp(lat_string, "obs", 3) == 0) { *plat_ref = b0; lat_stat = 1; } else lat_stat = sscanf(lat_string, "%f", plat_ref); if (seg_cmp(x0_string, "nom", 3) == 0) { *px0_out = nxout/2.; x0_stat = 1; } else x0_stat = sscanf(x0_string, "%f", px0_out); if (seg_cmp(y0_string, "nom", 3) == 0) { *py0_out = nyout/2.; y0_stat = 1; } else y0_stat = sscanf(y0_string, "%f", py0_out); if (lat_stat == 1 && lon_stat == 1 && x0_stat == 1 && y0_stat == 1) return 1; } } return EOF; } /* else parameters contained in program arguments. */ if (start == 0) return EOF; start = 0; if (seg_cmp(argv[3], "obs", 3) == 0) { *plon_ref = l0; lon_stat = 1; } else lon_stat = sscanf(argv[3], "%f", plon_ref); if (seg_cmp(argv[4], "obs", 3) == 0) { *plat_ref = b0; lat_stat = 1; } else lat_stat = sscanf(argv[4], "%f", plat_ref); if (seg_cmp(argv[5], "nom", 3) == 0) { *px0_out = nxout/2; x0_stat = 1; } else x0_stat = sscanf(argv[5], "%f", px0_out); if (seg_cmp(argv[6], "nom", 3) == 0) { *py0_out = nyout/2; y0_stat = 1; } else y0_stat = sscanf(argv[6], "%f", py0_out); if (sscanf(argv[7], "%s", out_file) == 1 && lon_stat == 1 && lat_stat == 1 && x0_stat == 1 && y0_stat == 1) return 1; return EOF; } /* Smallest power of 2 greater than or equal to the argument */ int lp2ge(int arg) { int i; for (i = 1; i < arg; i <<= 1); return i; } /* Function inset_data * This function rearranges image data originally formatted in a 2-D * array of size (nx1, ny1) into an array of size (nx2, ny2). */ int inset_data(float buff[], int nx1, int nx2, int ny1, int ny2) { int i, j; float temp, **dat1, **dat2, alpha, beta, arg, ds; if (nx2 <= nx1) return 0; dat1 = (float **)malloc(ny1*sizeof(float *)); dat2 = (float **)malloc(ny2*sizeof(float *)); for (i = 0; i < ny1; i++) dat1[i] = buff + nx1*i; for (i = 0; i < ny2; i++) dat2[i] = buff + nx2*i; for (i = ny1 - 1; i >= 1; i--) for (j = nx1 - 1; j >= 0; j--) { dat2[i][j] = dat1[i][j]; dat1[i][j] = 0.; } free(dat1); free(dat2); return 1; } /* Function expand_frame * This function redistributes datapoints so that the spacing is * mag pixels, where mag (magnification) is an external integer. * Pixels whose coordinates are other than (mag*i, mag*j) where * i and j are integers are set to 0. */ int expand_frame(float buff[]) { int i, j; float temp, **dat1, **dat2, alpha, beta, arg, ds; if (mag == 1) return 0; dat1 = (float **)malloc(nyframe*sizeof(float *)); dat2 = (float **)malloc(nymagframe*sizeof(float *)); for (i = 0; i < nyframe; i++) dat1[i] = buff + nxmagframe*i; for (i = 0; i < nyframe; i++) dat2[i] = buff + mag*nxmagframe*i; temp = dat1[0][0]; for (i = nyframe - 1; i >= 0; i--) for (j = nxframe - 1; j >= 0; j--) { dat2[i][mag*j] = dat1[i][j]; dat1[i][j] = 0.; } dat2[0][0] = temp; free(dat1); free(dat2); return 1; } /* Function read_data * This function reads the FITS file whose file ID is indicated by the * first argument and set the pointer indicated by the fourth argument * to the address an array that contains the array data. It is a * variation of the function "read_fits" in the file "do_fits.c". * The difference is that the amount for space allocated for the data * can be different (greater) than the amount of data read, as specified * by the third argument. The standard FITS dimensional and scaling * variables are read from the header and set accordingly. The FITS * file must be small enough that its contents can be fit into available * memory. */ int read_data(int fid, int n_dat, int n_buff, float **pbuff) { int i; *pbuff = (float *)malloc(n_buff*sizeof(float)); for (i = 0; i < n_buff; i++) (*pbuff)[i] = 0.; if (input_data(fid, *pbuff, n_dat) == EOF) { fprintf(stderr, "%s: Premature end of file encountered in read. ", dfprogname); fprintf(stderr, "Exiting.\n"); exit(0); } min_max(*pbuff, n_dat); return 1; } /* Function interp_data * This function Fourier interpolates data spread over an expanded * array, in both dimensions, filling in pixels whose coordinates are * other than (mag*i, mag*j), where i and j are integers. */ int datstep, magdatstep; int interp_data(float *buff, int magarg) { int i, j; /* First interpolate in the x-direction. */ datstep = 1; magdatstep = magarg; for (i = 0; i < nyframe; i++) f_interp(nxmagframe, buff + i*magarg*nxmagframe); /* Then interpolate in the y-direction. */ datstep = nxmagframe; magdatstep = magarg*nxmagframe; for (i = 0; i < nx; i++) f_interp(nymagframe, buff + i); return 1; } /* Function f_interp * This function Fourier interpolates data in a single dimension. */ int im_part_offset; int f_interp(int num_pts, float dat[]) { int i, j, num_grid, num_pairs; float theta, dwave_real, dwave_imag, wave_real, wave_imag; float temp, *ftrans_real, *ftrans_imag; num_grid = num_pts/mag; num_pairs = num_grid/2; ftrans_real = (float *)malloc(num_grid*sizeof(float)); im_part_offset = num_grid/2; ftrans_imag = ftrans_real + im_part_offset; rfft(num_grid, dat, ftrans_real, FORWARD); for (i = 1; i < mag; i++) { theta = pi/(double)(num_pairs*mag); dwave_real = sin(0.5*theta); dwave_real *= -2.*dwave_real; dwave_imag = -sin(theta); wave_real = 1.0; wave_imag = 0.0; for (j = 1; j < num_pairs; j++) { wave_real += dwave_real*(temp = wave_real) - dwave_imag*wave_imag; wave_imag += dwave_real*wave_imag + dwave_imag*temp; ftrans_real[j] = wave_real*(temp = ftrans_real[j]) - wave_imag*ftrans_imag[j]; ftrans_imag[j] = wave_real*ftrans_imag[j] + wave_imag*temp; } rfft(num_grid/2, ftrans_real, dat + i*datstep, INVERSE); } free(ftrans_real); return 1; } /* This function takes the Fourier transform of real data in a single * dimension whose order is a power of 2. */ #define STEP_REAL magdatstep #define STEP_COMPLEX 1 #define OFFSET im_part_offset int rfft(int n, float dat_in[], float dat_out[], int mode) { int i, j, numbits; int num_pix, num_pairs; float *in_real, *in_imag, *out_real, *out_imag; float even_real, even_imag, odd_real, odd_imag, theta, temp; float dwave_real, dwave_imag, wave_real, wave_imag; float recip_root_two, recip_two_root_two; recip_root_two = 1./sqrt(2.); recip_two_root_two = .5*recip_root_two; if ((numbits = count_bits(n)) < 0) return WRONG_NUM_COLS; if (mode == FORWARD) num_pix = n; else num_pix = 2*n; num_pairs = num_pix/2; in_real = (float *)malloc(num_pix*sizeof(float)); in_imag = in_real + num_pairs; out_real = (float *)malloc(num_pix*sizeof(float)); out_imag = out_real + num_pairs; if (mode == FORWARD) { for (j = 0; j < num_pairs; j++) { in_real[j] = dat_in[STEP_REAL* 2*j ]; in_imag[j] = dat_in[STEP_REAL*(2*j + 1)]; } if (fft(num_pairs, in_real, FORWARD) == EOF) return ROW_FAILURE; theta = pi/(double)num_pairs; dwave_real = sin(0.5*theta); dwave_real *= -2.*dwave_real; dwave_imag = sin(theta); dat_out[ 0] = in_real[0] + in_imag[0]; dat_out[OFFSET] = in_real[0] - in_imag[0]; wave_real = 1.0; wave_imag = 0.0; for (j = 1; j < num_pairs; j++) { even_real = recip_two_root_two*(in_real[j] + in_real[num_pairs - j]); even_imag = recip_two_root_two*(in_imag[j] - in_imag[num_pairs - j]); odd_real = recip_two_root_two*(in_imag[j] + in_imag[num_pairs - j]); odd_imag = -recip_two_root_two*(in_real[j] - in_real[num_pairs - j]); wave_real += dwave_real*(temp = wave_real) - dwave_imag*wave_imag; wave_imag += dwave_real*wave_imag + dwave_imag*temp; odd_real = wave_real*(temp = odd_real) - wave_imag*odd_imag; odd_imag = wave_real*odd_imag + wave_imag*temp; dat_out[ STEP_COMPLEX*j] = even_real + odd_real; dat_out[OFFSET + STEP_COMPLEX*j] = even_imag + odd_imag; } } else { for (j = 0; j < num_pairs; j++) { in_real[j] = dat_in[ STEP_COMPLEX*j]; in_imag[j] = dat_in[OFFSET + STEP_COMPLEX*j]; } theta = pi/(double)num_pairs; dwave_real = sin(0.5*theta); dwave_real *= -2.*dwave_real; dwave_imag = -sin(theta); out_real[0] = .5*(in_real[0] + in_imag[0]); out_imag[0] = .5*(in_real[0] - in_imag[0]); wave_real = 1.0; wave_imag = 0.0; for (j = 1; j < n; j++) { even_real = recip_root_two*(in_real[j] + in_real[num_pairs - j]); even_imag = recip_root_two*(in_imag[j] - in_imag[num_pairs - j]); odd_real = recip_root_two*(in_real[j] - in_real[num_pairs - j]); odd_imag = recip_root_two*(in_imag[j] + in_imag[num_pairs - j]); wave_real += dwave_real*(temp = wave_real) - dwave_imag*wave_imag; wave_imag += dwave_real*wave_imag + dwave_imag*temp; odd_real = wave_real*(temp = odd_real) - wave_imag*odd_imag; odd_imag = wave_real*odd_imag + wave_imag*temp; out_real[j] = even_real - odd_imag; out_imag[j] = even_imag + odd_real; } if (fft(n, out_real, INVERSE) == EOF) return ROW_FAILURE; for (j = 0; j < num_pairs; j++) { dat_out[STEP_REAL* 2*j ] = out_real[j]; dat_out[STEP_REAL*(2*j + 1)] = out_imag[j]; } } free(in_real); free(out_real); return 1; } int count_bits(int npts) { int i, num_bits, sign; for (num_bits = 0, i = npts, sign = 1; i != 1; i >>= 1) { if (i & 1 == 1) sign = -1; num_bits++; } return sign*num_bits; } /* Function to fast_Fourier transform the array, dat. */ /* This function inputs an array of float pairs (complex) */ /* of length npts (int). The Fourier transform of dat replaces */ /* the original contents of dat. */ int fft(int npts, float dat[], int mode) { unsigned int i, j, k_e, k_o, seg_length; double theta, wave_real, wave_imag, dwave_real, dwave_imag; float factor, tmp_real, tmp_imag, *dat_real, *dat_imag; /* If npts is not a power of 2 the value -1, indicating an error */ /* condition, is returned and no Fourier transform is attempted. */ dat_real = dat; dat_imag = dat + npts; bit_rev(npts, dat); for (seg_length = 1; seg_length < npts; seg_length <<= 1) { theta = pi/(double)seg_length; dwave_real = sin(0.5*theta); dwave_real *= -2.*dwave_real; if (mode >= 0) dwave_imag = sin(theta); else dwave_imag = -sin(theta); wave_real = 1.0; wave_imag = 0.0; for (j = 0; j < seg_length; j++) { for (k_e = j; k_e < npts; k_e += seg_length << 1) { k_o = k_e + seg_length; tmp_real = wave_real*dat_real[k_o] - wave_imag*dat_imag[k_o]; tmp_imag = wave_imag*dat_real[k_o] + wave_real*dat_imag[k_o]; dat_real[k_o] = dat_real[k_e] - tmp_real; dat_imag[k_o] = dat_imag[k_e] - tmp_imag; dat_real[k_e] += tmp_real; dat_imag[k_e] += tmp_imag; } wave_real += (tmp_real = wave_real)*dwave_real - wave_imag*dwave_imag; wave_imag += tmp_real*dwave_imag + wave_imag*dwave_real; } } factor = 1./(float)sqrt((double)npts); for (i = 0; i < npts; i++) { dat_real[i] *= factor; dat_imag[i] *= factor; } return 1; } /* Function bit_rev * This function performs bit-reversal remapping of data represented by * its second argument. */ int bit_rev(int npts, float dat[]) { int i, j, arg_i, arg_j, k, kstart; float *dat_real, *dat_imag, tmp; dat_real = dat; dat_imag = dat + npts; for (i = (j = 0), kstart = npts >> 1; i < npts; i++) { if (j > i) { arg_i = i; arg_j = j; tmp = dat_real[arg_i]; dat_real[arg_i] = dat_real[arg_j]; dat_real[arg_j] = tmp; tmp = dat_imag[arg_i]; dat_imag[arg_i] = dat_imag[arg_j]; dat_imag[arg_j] = tmp; } for (k = kstart; k > 0 && j >= k; k >>= 1) j -= k; j += k; } return 1; } /* Function los_project * This function uses cubic interpolation to remap data by line-of-sight * projection. */ int los_project(float lon_ref, float lat_ref, float buff_in[], float buff_out[]) { int i, j; float x, y, z, b0rad, **dat_in, **dat_out; float temp; float lon_rad, lat_rad, p_angle, psi; float cos_b0, sin_b0, cos_p, sin_p, cos_l, sin_l, cos_b, sin_b, cos_2psi, sin_2psi; float eta, new_eta, xi, new_xi, zeta, zeta2; float rot_xx, rot_xy, rot_xz, rot_yx, rot_yy, rot_yz, rot_zx, rot_zy, rot_zz; float squeeze_xx, squeeze_xy, squeeze_yx, squeeze_yy; float eta2, xi2, rho, phi; float sin_alpha, cos_alpha, tan_alpha, tan_ratio, parallax; float sin_alpha_max, cos_alpha_max, tan_alpha_max, sin2_diff; b0rad = b0*rad_fact; p_angle = pangle*rad_fact; lon_rad = (lon_ref - l0)*rad_fact; lat_rad = lat_ref*rad_fact; psi = maj_ax_ang*rad_fact; sin_alpha_max = sin(sun_rad_radian); cos_alpha_max = cos(sun_rad_radian); tan_alpha_max = sin_alpha_max/cos_alpha_max; cos_b0 = cos(-b0rad); sin_b0 = sin(-b0rad); cos_p = cos(-p_angle); sin_p = sin(-p_angle); cos_l = cos(-lon_rad); sin_l = sin(-lon_rad); cos_b = cos( lat_rad); sin_b = sin( lat_rad); cos_2psi = cos( 2.*psi); sin_2psi = sin( 2.*psi); dat_in = (float **)malloc(nymagframe*sizeof(float *)); dat_out = (float **)malloc(nyout*sizeof(float *)); for (i = 0; i < nymagframe; i++) dat_in[i] = buff_in + nxmagframe*i; for (i = 0; i < nyout; i++) dat_out[i] = buff_out + nxout*i; rot_xx = cos_p*cos_l + sin_p* sin_b0* sin_l ; rot_xy = cos_p*sin_l*sin_b + sin_p*(cos_b0*cos_b - sin_b0*sin_b*cos_l); rot_xz = -cos_p*sin_l*cos_b + sin_p*(cos_b0*sin_b + sin_b0*cos_b*cos_l); rot_yx = -sin_p*cos_l + cos_p* sin_b0* sin_l ; rot_yy = -sin_p*sin_l*sin_b + cos_p*(cos_b0*cos_b - sin_b0*sin_b*cos_l); rot_yz = sin_p*sin_l*cos_b + cos_p*(cos_b0*sin_b + sin_b0*cos_b*cos_l); rot_zx = cos_b0* sin_l ; rot_zy = -sin_b0*cos_b - cos_b0*sin_b*cos_l ; rot_zz = -sin_b0*sin_b + cos_b0*cos_b*cos_l ; squeeze_xx = .5*parity*((sun_rad_maj + sun_rad_min) - (sun_rad_maj - sun_rad_min)*cos_2psi); squeeze_xy = .5*parity* (sun_rad_maj - sun_rad_min)*sin_2psi; squeeze_yx = .5 * (sun_rad_maj - sun_rad_min)*sin_2psi; squeeze_yy = .5 *((sun_rad_maj + sun_rad_min) + (sun_rad_maj - sun_rad_min)*cos_2psi); /* In the following loop, * xi is the distance in the y-direction from frame center; * eta is the distance in the x-direction from frame center; * xi and eta in the output frame correspond to x and y in the * input frame. */ for (i = 0; i < nyout; i++) /* eta index */ { eta = dyout*((float)i - y_center_out); eta2 = eta*eta; for (j = 0; j < nxout; j++) /* xi index */ { xi = dxout*((float)j - x_center_out); xi2 = xi*xi; if (eta == 0. && xi == 0.) phi = 0.; else phi = atan2(eta, xi); tan_ratio = sqrt(xi2 + eta2); tan_alpha = tan_alpha_max*tan_ratio; cos_alpha = 1./sqrt(1. + tan_alpha*tan_alpha); sin_alpha = tan_alpha*cos_alpha; sin2_diff = sin_alpha_max*sin_alpha_max - sin_alpha*sin_alpha; if (sin2_diff > 0.) { rho = (cos_alpha - sqrt(sin2_diff)) * (sin_alpha/sin_alpha_max); zeta = sqrt(1. - rho*rho); new_xi = rho*cos(phi); new_eta = rho*sin(phi); z = rot_zx*new_xi + rot_zy*new_eta + rot_zz*zeta; if (z > sin_alpha_max) { x = rot_xx*new_xi + rot_xy*new_eta + rot_xz*zeta; y = rot_yx*new_xi + rot_yy*new_eta + rot_yz*zeta; parallax = sqrt(1. - sin_alpha_max*sin_alpha_max) / (1. - sin_alpha_max*z); x *= parallax; y *= parallax; temp = x; x = x_center + squeeze_xx*temp + squeeze_xy*y; y = y_center + squeeze_yx*temp + squeeze_yy*y; dat_out[i][j] = bicubic_interp(x, y, dat_in); } else dat_out[i][j] = DARK_SIDE_VALUE; } else dat_out[i][j] = SPACE_VALUE; } } free(dat_in); free(dat_out); return 1; } /* Function perp_project * This function uses cubic interpolation to remap data by perpendicular * projection. */ int perp_project(float lon_ref, float lat_ref, float buff_in[], float buff_out[]) { int i, j; float x, y, z, b0rad, **dat_in, **dat_out; float temp; float lon_rad, lat_rad, p_angle, psi; float cos_b0, sin_b0, cos_p, sin_p, cos_l, sin_l, cos_b, sin_b, cos_2psi, sin_2psi; float eta, new_eta, xi, new_xi, zeta, zeta2; float rot_xx, rot_xy, rot_xz, rot_yx, rot_yy, rot_yz, rot_zx, rot_zy, rot_zz; float squeeze_xx, squeeze_xy, squeeze_yx, squeeze_yy; float eta2, xi2, rho, rho2, phi; float sin_alpha, cos_alpha, tan_alpha, parallax; float sin_alpha_max, cos_alpha_max, tan_alpha_max; b0rad = b0*rad_fact; p_angle = pangle*rad_fact; lon_rad = (lon_ref - l0)*rad_fact; lat_rad = lat_ref*rad_fact; psi = maj_ax_ang*rad_fact; sin_alpha_max = sin(sun_rad_radian); cos_b0 = cos(-b0rad); sin_b0 = sin(-b0rad); cos_p = cos(-p_angle); sin_p = sin(-p_angle); cos_l = cos(-lon_rad); sin_l = sin(-lon_rad); cos_b = cos( lat_rad); sin_b = sin( lat_rad); cos_2psi = cos( 2.*psi); sin_2psi = sin( 2.*psi); dat_in = (float **)malloc(nymagframe*sizeof(float *)); dat_out = (float **)malloc(nyout*sizeof(float *)); for (i = 0; i < nymagframe; i++) dat_in[i] = buff_in + nxmagframe*i; for (i = 0; i < nyout; i++) dat_out[i] = buff_out + nxout*i; rot_xx = cos_p*cos_l + sin_p* sin_b0* sin_l ; rot_xy = cos_p*sin_l*sin_b + sin_p*(cos_b0*cos_b - sin_b0*sin_b*cos_l); rot_xz = -cos_p*sin_l*cos_b + sin_p*(cos_b0*sin_b + sin_b0*cos_b*cos_l); rot_yx = -sin_p*cos_l + cos_p* sin_b0* sin_l ; rot_yy = -sin_p*sin_l*sin_b + cos_p*(cos_b0*cos_b - sin_b0*sin_b*cos_l); rot_yz = sin_p*sin_l*cos_b + cos_p*(cos_b0*sin_b + sin_b0*cos_b*cos_l); rot_zx = cos_b0* sin_l ; rot_zy = -sin_b0*cos_b - cos_b0*sin_b*cos_l ; rot_zz = -sin_b0*sin_b + cos_b0*cos_b*cos_l ; squeeze_xx = .5*parity*((sun_rad_maj + sun_rad_min) - (sun_rad_maj - sun_rad_min)*cos_2psi); squeeze_xy = .5*parity* (sun_rad_maj - sun_rad_min)*sin_2psi; squeeze_yx = .5 * (sun_rad_maj - sun_rad_min)*sin_2psi; squeeze_yy = .5 *((sun_rad_maj + sun_rad_min) + (sun_rad_maj - sun_rad_min)*cos_2psi); /* In the following loop, * xi is the distance in the y-direction from frame center; * eta is the distance in the x-direction from frame center; * xi and eta in the output frame correspond to x and y in the * input frame. */ for (i = 0; i < nyout; i++) /* eta index */ { eta = dyout*((float)i - y_center_out); eta2 = eta*eta; for (j = 0; j < nxout; j++) /* xi index */ { xi = dxout*((float)j - x_center_out); xi2 = xi*xi; if (eta == 0. && xi == 0.) phi = 0.; else phi = atan2(eta, xi); rho2 = xi2 + eta2; rho = sqrt(rho2); zeta2 = 1. - rho2; if (zeta2 >= 0.) zeta = sqrt(zeta2); else zeta = -1.; new_xi = rho*cos(phi); new_eta = rho*sin(phi); if (zeta > 0.) { z = rot_zx*new_xi + rot_zy*new_eta + rot_zz*zeta; if (z > sin_alpha_max) { x = rot_xx*new_xi + rot_xy*new_eta + rot_xz*zeta; y = rot_yx*new_xi + rot_yy*new_eta + rot_yz*zeta; parallax = sqrt(1. - sin_alpha_max*sin_alpha_max) / (1. - sin_alpha_max*z); x *= parallax; y *= parallax; temp = x; x = x_center + squeeze_xx*temp + squeeze_xy*y; y = y_center + squeeze_yx*temp + squeeze_yy*y; dat_out[i][j] = bicubic_interp(x, y, dat_in); } else dat_out[i][j] = DARK_SIDE_VALUE; } else dat_out[i][j] = SPACE_VALUE; } } free(dat_in); free(dat_out); return 1; } /* Function postel * This function uses cubic interpolation to remap data by Postel * projection. */ int postel(float lon_ref, float lat_ref, float buff_in[], float buff_out[]) { int i, j; float x, y, z, b0rad, **dat_in, **dat_out; float temp; float lon_rad, lat_rad, p_angle, psi; float cos_b0, sin_b0, cos_p, sin_p, cos_l, sin_l, cos_b, sin_b, cos_2psi, sin_2psi; float eta, new_eta, xi, new_xi, zeta, zeta2; float rot_xx, rot_xy, rot_xz, rot_yx, rot_yy, rot_yz, rot_zx, rot_zy, rot_zz; float squeeze_xx, squeeze_xy, squeeze_yx, squeeze_yy; float eta2, xi2, rho, rho2, phi; float sin_alpha_max, parallax; b0rad = b0*rad_fact; p_angle = pangle*rad_fact; lon_rad = (lon_ref - l0)*rad_fact; lat_rad = lat_ref*rad_fact; psi = maj_ax_ang*rad_fact; sin_alpha_max = sin(sun_rad_radian); cos_b0 = cos(-b0rad); sin_b0 = sin(-b0rad); cos_p = cos(-p_angle); sin_p = sin(-p_angle); cos_l = cos(-lon_rad); sin_l = sin(-lon_rad); cos_b = cos( lat_rad); sin_b = sin( lat_rad); cos_2psi = cos( 2.*psi); sin_2psi = sin( 2.*psi); dat_in = (float **)malloc(nymagframe*sizeof(float *)); dat_out = (float **)malloc(nyout*sizeof(float *)); for (i = 0; i < nymagframe; i++) dat_in[i] = buff_in + nxmagframe*i; for (i = 0; i < nyout; i++) dat_out[i] = buff_out + nxout*i; rot_xx = cos_p*cos_l + sin_p* sin_b0* sin_l ; rot_xy = cos_p*sin_l*sin_b + sin_p*(cos_b0*cos_b - sin_b0*sin_b*cos_l); rot_xz = -cos_p*sin_l*cos_b + sin_p*(cos_b0*sin_b + sin_b0*cos_b*cos_l); rot_yx = -sin_p*cos_l + cos_p* sin_b0* sin_l ; rot_yy = -sin_p*sin_l*sin_b + cos_p*(cos_b0*cos_b - sin_b0*sin_b*cos_l); rot_yz = sin_p*sin_l*cos_b + cos_p*(cos_b0*sin_b + sin_b0*cos_b*cos_l); rot_zx = cos_b0* sin_l ; rot_zy = -sin_b0*cos_b - cos_b0*sin_b*cos_l ; rot_zz = -sin_b0*sin_b + cos_b0*cos_b*cos_l ; squeeze_xx = .5*parity*((sun_rad_maj + sun_rad_min) - (sun_rad_maj - sun_rad_min)*cos_2psi); squeeze_xy = .5*parity* (sun_rad_maj - sun_rad_min)*sin_2psi; squeeze_yx = .5 * (sun_rad_maj - sun_rad_min)*sin_2psi; squeeze_yy = .5 *((sun_rad_maj + sun_rad_min) + (sun_rad_maj - sun_rad_min)*cos_2psi); /* In the following loop, * xi is the distance in the y-direction from frame center; * eta is the distance in the x-direction from frame center; * xi and eta in the output frame correspond to x and y in the * input frame. */ for (i = 0; i < nyout; i++) /* eta index */ { eta = dyout*((float)i - y_center_out); eta2 = eta*eta; for (j = 0; j < nxout; j++) /* xi index */ { xi = dxout*((float)j - x_center_out); xi2 = xi*xi; if (eta == 0. && xi == 0.) phi = 0.; else phi = atan2(eta, xi); rho2 = xi2 + eta2; rho = sqrt(rho2); zeta = cos(rho); rho = sin(rho); new_xi = rho*cos(phi); new_eta = rho*sin(phi); z = rot_zx*new_xi + rot_zy*new_eta + rot_zz*zeta; if (z > sin_alpha_max) { x = rot_xx*new_xi + rot_xy*new_eta + rot_xz*zeta; y = rot_yx*new_xi + rot_yy*new_eta + rot_yz*zeta; parallax = sqrt(1. - sin_alpha_max*sin_alpha_max) / (1. - sin_alpha_max*z); x *= parallax; y *= parallax; temp = x; x = x_center + squeeze_xx*temp + squeeze_xy*y; y = y_center + squeeze_yx*temp + squeeze_yy*y; dat_out[i][j] = bicubic_interp(x, y, dat_in); } else dat_out[i][j] = DARK_SIDE_VALUE; } } free(dat_in); free(dat_out); return 1; } /* Function bicubic_interp * This function performs cubic interpolation in two dimensions. */ float bicubic_interp(float x, float y, float **dat_in) { int k, i_0, i_1, i_2, i_3, j_0, j_1, j_2, j_3; float alpha_x, alpha_y, interp_grid[4]; i_1 = gile(y); i_0 = i_1 - 1; i_2 = i_1 + 1; i_3 = i_2 + 1; j_1 = gile(x); j_0 = j_1 - 1; j_2 = j_1 + 1; j_3 = j_2 + 1; if (i_1 < 0 || j_1 < 0 || i_2 > ny - 1 || j_2 > nx - 1) return 0.; if (i_3 > ny - 1) { i_0--; i_1--; i_2--; i_3--; } if (j_3 > nx - 1) { j_0--; j_1--; j_2--; j_3--; } if (i_0 < 0) { i_0++; i_1++; i_2++; i_3++; } if (j_0 < 0) { j_0++; j_1++; j_2++; j_3++; } alpha_x = x - (float)j_1; alpha_y = y - (float)i_1; for (k = 0; k <= 3; k++) interp_grid[k] = cu_interp(dat_in[i_0 + k] + j_0, alpha_x); return cu_interp(interp_grid, alpha_y); } /* Function cu_interp * This function performs cubic interpolation in one dimension. */ float cu_interp(float ref_grid[], float x) { int i, j; float sum, xpwr; static float coeff_buff[] = { 0., 1., 0., 0., -1./3., -1./2., 1., -1./6., 1./2., -1., 1./2., 0., -1./6., 1./2., -1./2., 1./6. }; static float *coeff[] = { coeff_buff, coeff_buff + 4, coeff_buff + 8, coeff_buff + 12 }; for (i = 0, sum = 0., xpwr = 1.; i <= 3; i++, xpwr *= x) for (j = 0; j <= 3; j++) sum += ref_grid[j]*coeff[i][j]*xpwr; return sum; } /* Function gile * This function returnes the greatest integer less than or equal to * its floating-point argument. */ int gile(float x) { int i; i = (int)x; if (x >= 0.) return i; if ((float)i != x) return i - 1; return i; } /* Function pad_image * This function cyclically interpolates pixels outside of the limb along * horizontal lines connecting the limbs where horizontal lines intersect * the solar image. Outside of the band in which horizontal lines connect * to the solar limb, the interpolation is done vertically between the * bottom and top of the band. The purpose is to smooth out the limb * transition to minimize ringing of the discontinuity in response to the * Fourier interpolation. */ int pad_image(float data[], float xc, float yc, float rad) { int i, i1, i2, j, j1, j2; float thresh2; float xdist, xdist2, ydist, ydist2, *dat, *datref1, *datref2; float alpha, beta; thresh2 = rad*rad; i1 = (int)(yc - rad + .5); i2 = (int)(yc + rad + .5); if (i1 < 0) i1 = 0; if (i2 < 0) i2 = 0; if (i1 > ny - 1) i1 = ny - 1; if (i2 > ny - 1) i2 = ny - 1; for (i = i1; i <= i2; i++) { dat = data + i*nx; ydist = (float)i - yc; ydist2 = ydist*ydist; xdist2 = thresh2 - ydist2; if (xdist2 < 0.) xdist = 0.; else xdist = sqrt(xdist2); j1 = (int)(xc - xdist + .5); j2 = (int)(xc + xdist + .5); if (j1 < 0) j1 = 0; if (j2 < 0) j2 = 0; if (j1 > nx - 1) j1 = nx - 1; if (j2 > nx - 1) j2 = nx - 1; for (j = 0; j < j1; j++) { beta = (float)(j1 - j)/(float)(j1 + nx - j2); alpha = 1. - beta; dat[j] = alpha*dat[j1] + beta*dat[j2]; } for (j = j2 + 1; j < nx; j++) { beta = (float)(j1 - j + nx)/(float)(j1 + nx - j2); alpha = 1. - beta; dat[j] = alpha*dat[j1] + beta*dat[j2]; } } for (i = 0; i < i1; i++) { beta = (float)(i1 - i)/(float)(i1 + ny - i2); dat = data + i *nx; datref1 = data + i1*nx; datref2 = data + i2*nx; for (j = 0; j < nx; j++) dat[j] = alpha*datref1[j] + beta*datref2[j]; } for (i = i2 + 1; i < ny; i++) { beta = (float)(i1 - i + ny)/(float)(i1 + ny - i2); dat = data + i *nx; datref1 = data + i1*nx; datref2 = data + i2*nx; for (j = 0; j < nx; j++) dat[j] = alpha*datref1[j] + beta*datref2[j]; } return 1; } /* This function truncates the value of the image to a constant value * some distance inside of the actual limb. */ int truncate_limb(float data_in[], float xc, float yc, float trunc_rad, float value) { int i, j; float *dat; float xdist, ydist, ydist2, rad; for (i = 0; i < ny; i++) { dat = data_in + nx*i; ydist = (float)i - yc; ydist2 = ydist*ydist; for (j = 0; j < nx; j++) { xdist = (float)j - xc; rad = sqrt(xdist*xdist + ydist2); if (rad > trunc_rad) dat[j] = value; } } return 1; } /* Function meas_rot * This function fits the data inside of a limb threshold to a planar * function of position and returns the fits parameters to addresses * to which the last five arguments point. */ int meas_rot(float data[], float *prefx, float *prefy, float *pmean, float *pgradx, float *pgrady) { int i, j, k, index; float refx, refy, mean_xx, mean_yy, mean_xy; float mean_dat, mean_xdat, mean_ydat; float x, x2, y, y2, r, r2, factor, rthresh, rthresh2; rthresh = sun_rad_min*(1. - ROT_STATIST_MARGIN); rthresh2 = rthresh*rthresh; refx = 0.; refy = 0.; mean_dat = 0.; for (i = 0, k = 0; i < ny; i++) { y = i - y_center; y2 = y*y; for (j = 0; j < nx; j++) { x = j - x_center; x2 = x*x; r2 = x2 + y2; index = nx*i + j; if (r2 < rthresh2) { refx += x; refy += y; mean_dat += data[index]; k++; } } } if (k == 0) { *prefx = 0.; *prefy = 0.; *pmean = 0.; *pgradx = 0.; *pgrady = 0.; return EOF; } factor = 1./k; refx *= factor; refy *= factor; mean_dat *= factor; mean_xx = 0.; mean_xy = 0.; mean_yy = 0.; mean_xdat = 0.; mean_ydat = 0.; for (i = 0; i < ny; i++) { y = i - y_center; y2 = y*y; y -= refy; for (j = 0; j < nx; j++, k++) { x = j - x_center; x2 = x*x; r2 = x2 + y2; x -= refx; index = nx*i + j; if (r2 < rthresh2) { mean_xx += x*x; mean_xy += x*y; mean_yy += y*y; mean_xdat += x*data[index]; mean_ydat += y*data[index]; } } } mean_xx *= factor; mean_xy *= factor; mean_yy *= factor; mean_xdat *= factor; mean_ydat *= factor; if ((factor = mean_xx*mean_yy - mean_xy*mean_xy) <= 0.) { *prefx = 0.; *prefy = 0.; *pmean = 0.; *pgradx = 0.; *pgrady = 0.; return EOF; } factor = 1./factor; *prefx = refx; *prefy = refy; *pmean = mean_dat; *pgradx = (mean_xdat*mean_yy - mean_ydat*mean_xy)*factor; *pgrady = (mean_ydat*mean_xx - mean_xdat*mean_xy)*factor; return 1; } /* Function remove_rot * This function inputs arguments that characterize solar rotation and * orbital motion and subtracts these from the data. */ float remove_rot(float data[], float refx, float refy, float mean, float gradx, float grady) { int i, j, index, n_stat; float x, x2, y, y2, r, r2; float rthresh, rthresh2; float rthresh_stat, rthresh_stat2; float sum_sq, rms; rthresh = sun_rad_maj; rthresh2 = rthresh*rthresh; rthresh_stat = sun_rad_min*(1. - ROT_STATIST_MARGIN); rthresh_stat2 = rthresh_stat*rthresh_stat; for (i = 0, sum_sq = 0., n_stat = 0; i < ny; i++) { y = (float)i - y_center; y2 = y*y; y -= refy; for (j = 0; j < nx; j++) { x = (float)j - x_center; x2 = x*x; x -= refx; r2 = x2 + y2; index = nx*i + j; if (r2 < rthresh2) data[index] -= mean + x*gradx + y*grady; else data[index] = 0.; if (r2 < rthresh_stat2) { sum_sq += data[index]*data[index]; n_stat++; } } } sum_sq /= n_stat; rms = sqrt(sum_sq); return rms; } /* This function interprets the value of the ITIME keyword in the * MDI image and returns it as a floating-point number of seconds * to the address to which its second argument points. */ int getdt(char line[], float *pdtime) { int i; float unit; unit = 1.; for (i = 0; i < LINE_LENGTH; i++) { if (line[i] == '\0') break; if ((line[i] >= '0' && line[i] <= '9') || line[i] == '-' || line[i] == '.') ; else { if (line[i] == 's') { unit = 1.0; line[i] = '\0'; break; } else if (line[i] == 'm') { unit = 60.0; line[i] = '\0'; break; } else if (line[i] == 'h') { unit = 3600.0; line[i] = '\0'; break; } else line[i] = ' '; } } sscanf(line, "%e", pdtime); *pdtime *= unit; return 1; } /* This function interprets the value of the T_OBS keyword in the * MDI image and returns it as an an integer number of seconds since * 1996.0. */ int gettime(char line[], int *ptime) { int i, year, month, day, hr, min, sec, num_days; for (i = 0; i < LINE_LENGTH; i++) { if (line[i] == '\0') break; if (line[i] >= '0' && line[i] <= '9') ; else line[i] = ' '; } if (sscanf(line, "%d %d %d %d %d %d", &year, &month, &day, &hr, &min, &sec) != 6) return(EOF); num_days = date_to_jday(year, month, day); *ptime = sec + 60*(min + 60*(hr + 24*num_days)); return 1; } /* #define REF_DAY 1721060 */ /* Use this for Julian day */ #define REF_DAY (-729023) /* Days since 1996.0 */ #define NUM_MONTHS 12 /* Function date_to_jday * This function inputs the year, month, and day as arguments and computes * the number of days since 1996.0. */ int date_to_jday(int year, int month, int dom) { int nleaps, dfsoytsom; static int nomcumdays[] = { 0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334, 365}; static int leapcumdays[] = { 0, 31, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335, 366}; nleaps = (year - 1)/4 - year/100 + year/400; if (isleap(year) == 1) dfsoytsom = leapcumdays[month - 1]; /* days from start of year */ else /* to start of month */ dfsoytsom = nomcumdays[month - 1]; return 365*year + nleaps + dfsoytsom + dom + REF_DAY; } /* This function inputs the the number of days since 1996.0 and writes * the alpha-numeric date to the string to which its second argument * points. */ static char *moid[] = { "Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec" }; static int nomcumdays[] = { 0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334, 365}; static int leapcumdays[] = { 0, 31, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335, 366}; int jday_to_date(int day, char alpha_date[], int alpha_set) { int n1, n2, n3, n4, m1, m2, m3, m4; int i, year, month, dom, dday; dday = day - REF_DAY; if (dday < 366) { fprintf(stderr, "Sorry, I can only go back as far as %d.\n", REF_DAY + 366); exit(0); } n1 = dday/146097; /* 4-century number */ m1 = dday%146097; if (m1 < 36525) /* century 0 */ { n2 = 0; /* century number */ m2 = m1; n3 = m2/1461; /* 4-yr number */ m3 = m2%1461; if (m3 < 366) { n4 = m3/366; /* year number */ m4 = m3 + 1; } else { n4 = (m3 - 1)/365; m4 = (m3 - 1)%365 + 1; } } else /* centuries 1--3 (skip 1st leap-yr) */ { n2 = (m1 - 1)/36524; /* century number */ m2 = (m1 - 1)%36524; if (m2 < 1460) /* no leap year */ { n3 = 0; /* 4-yr number */ m3 = m2; n4 = m3/365; m4 = m3%365 + 1; } else /* leap year */ { n3 = (m2 + 1)/1461; /* 4-yr number */ m3 = (m2 + 1)%1461; if (m3 < 366) { n4 = m3/366; /* year number */ m4 = m3 + 1; } else { n4 = (m3 - 1)/365; m4 = (m3 - 1)%365 + 1; } } } year = 400*n1 + 100*n2 + 4*n3 + n4; if (isleap(year) == 1) { for (i = 0; i <= NUM_MONTHS; i++) if (leapcumdays[i] >= m4) break; month = --i; dom = m4 - leapcumdays[month]; } else { for (i = 0; i <= NUM_MONTHS; i++) if (nomcumdays[i] >= m4) break; month = --i; dom = m4 - nomcumdays[month]; } if (alpha_set == 1) sprintf(alpha_date, "%d %s %02d", year, moid[month], dom); else sprintf(alpha_date, "%d %02d %02d", year, month + 1, dom); return 1; } /* This function inputs the integer year as an argument. It returns * a value of 1 to identify a leap year and a value of 0 otherwise. */ int isleap(int year) { if (year%4 == 0 && (year%100 != 0 || year%400 == 0)) return(1); return 0; }